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ABSTRACT 


This work presents a comprehensive computer simulation 
for the hydrolytic polymerization of £ - caprolactam without 
water removal in a continuous tubular reactor. This incorporates 
the existence of a laminar ( Hagen-Poiseuille) velocity profile 
and radial thermal diffusion. Solution of the balance equations 
are carried out using two types of finite difference technique® 
as well as the method of orthogonal collocation coupled with 
GEAR'S method for ODE s. The simulation results are obtained 
for isothermal/ adiabatic and non-adiabatic operations for both 
plug and laminar flow velocity profiles. Some of the results 
were compared with available results for simpler cases. The 
influence of various operational variables on the reactor 
performance was studied. A comparison is made of the different 
methods used and their limitations discussed, ihe first method 
of finite difference (method I) was found to be superior to 
the second one, though both of them gave the same results for 
higher- number of grid points. The orthogonal collocation 
technique gave good results except for the laminar flow 
adiabatic and non-adiabatic cases where the temperature 
gradients near the wall were s^eep. An improvenent is 
suggested by including the finite element technique. 



CHAPTER 1 


INTRODUCTICa^ 


The hydrolytic polymerization of t -caprolactam is a very 

important commercial process and has drawn the attention of 

various researchers in recent years . in fact several excellent 
1-5 

reviews have appeared on this subject, emphasizing various 
aspects of the polymerization process. 

The worlc on nylon-6 polymerization falls broadly into two 

categories. The emphasis in the early work was primarily on the 

determination of the polymerization mechanism and of the rate and 

equilibritun constants for the various reactions. Ihese have been 

reviewed by Reimschuessel^. More recently the emphasis has shifted 

to the study of physical processes like heat transfer, mass transfer 

etc., associated with chemical reaction in large scale reactors, 

and to the actual modeling and optimization of industrial reactors. 

2 3 4 

Tai and Tagawa , and Kumar and Gupta ' have reviewed these 
aspects*. 

The reaction mechanian comprises of three major reactions, 
namely, the ring opening of caprolactam by water to form amino- 
caproic acid, and chain and step growth reactions, in addition, 
there are several important side reactions. Among these is the 
formation of cyclic oligomers, desamination and peroxidation of 
caprolactam^. The most important side reactions are those' 
associated with cyclic oligomers, since their presence in the 
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product causes problems in its p 2 x>cessing (e.g. , in spinning 
and molding) . Tbe kinetic scheme considered in this work 
is shown in Table I. Ihis includes the main reactions and 
the reactions associated with the cyclic dimer. The other 
cyclization reactions are omitted since their rate constants 
are not yet available. Moreover, it is well known ' that 
the formation of the cyclic dimer predominates, and can be 
used as a first-order approximation of the total cyclics 
present. 

The reactions are known to be catalyzed by the carboxyl 
end groups present in the reaction mass. The apparent rate 
constants are of the form 

= k° + [ -COOH] (1) 

with Arrhenius forms being used for k? and . All the 

reactions in Table I are reversible in nature end the 

temperature dependency of the equilibrium constants are 

given by standard thermodynamic relations. The rate and 

equilibrium constants for the reactions of Table I are given 

in Table II, These are based on results from a series of 

6 

esqseriments carried out by Tai et al., using a nonlinear 
regression analysis. 

Various kinds of reactors, e.g, , batch, tubular, 
con tinuous-flow- stirred- tank reactors (CSTRs) and combinations 
thereof, are employed in industry to manufacture nylon-6, 
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Table I 


Kinetic Scheme For Nylon-6 Polymerization 


1 . Ring Opening : 


+ W 


= k^/K^ 


2. Poly condensation: 

k„ 


p + P 
n m 


^2 


^n-HT, + w , n, p = 1,2,3 


3, Polyaddition: 


P_ + C, 


n • n ^^5 ’^n+l ' 


4. Ring Opening of Cyclic Dimer j 

k. 


C2 + W ~ 


"4 




k^ = k^/K^ 


5, Polyaddition of Cyclic Dimer: 

k. 


P_ + c 


n ^2 ^ 


^5 ^5 


^n+2 ' ” 


Contd 




Table I ( Continued) 


H-N- ( CH. 

I : 


H-N- ( CH. 


H 

H ^ N - 



( e - caprolactam) 



H 

N - (CH2)5 


0 

II 

- C ( cyclic dimer) 


J 


( CH2 ) ^ 



OH 


{ Polymer : Ny lon-6 ) 



Ecrui librium Constants for Nylon~6 Polymerization 


5 
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Computer simulations based on the mathematical modeling 

of the polymerization processes in these reactors or their 

combinations offer information which is of paramount 

importance for quality control, process control and 

operational optimization of existing plants as well as in 

the design of new plants. Several studies have been 

7-9 

reported on the simulation of both ideal reactors and 
some common industrial reactors^*^^^. 

Ihe present work focusses on the simulation of 
continuous tubular reactors with radial variations of 
temperature and concentrations accounted for. In most of 
the studies on tubular reactors reported till now, radial 
variations of temperature (and hence, of the concentrations) 
have not been considered, Ihis was probably to keep the 
analysis simple. The only study wherein such radial profiles 

7 

are indeed computed is the one by Tai et al. . '^hese 
workers have used the two dimensional (in the radial and 
axial directions) finite difference computational technique 
to solve their balance equations. However, they have also 
presented results on several other types of reactors and 
reactor combinations, and so their work on tubular reactors 
is slightly less comprehensive than what one would like to 
have. Moreover, our initial attempts on the optimization of 
such reactors has revealed the need for faster and more 
efficient computational techniques to integrate the modeling 
equations for such reactors, in order to cut down on the 



total computational costs. This work was undertaken to 
satisfy these needs, and presents simulation results using 
two types of finite difference algorithms, as well as th e 
orthogonal collocation technique. Results are generated 
for several conditions, e.g,, for both plug and laminar 
flow velocity profiles in the tubular reactor and for both 
adiabatic and non-adiabatic operations, A detailed 
comparison is made of the results obtained with various 
computational techniques and model assumptions. 



CHAPTER 2 


FORMULATICaa 


The reaction mechanism considered here is given in 
Table I with the rate and equilibrium constants as shown 
in Table II. The rates of formation of various ^ecies 
in the reaction mixture and those of the first two moments, 

Xq and are listed in Table III along with the rate of 

heat generation. In general, the k moment is defined as: 

= 0 / 1^2 ( 2 ) 

n=l " 

The widely accepted use of moments instead of 

individual concentrations of polymer molecules of various 

chain lengths helps reduce the number of equations ( for 

the individual species) to be solved simultaneously, and 

3 17—19 

gives almost the same results * . At the same time, 

information about the first two or three moments are usually 
sufficient to know the nature of the molecular weight 
distribution in the product polymer. For example, the 
number average chain length (Fjj) obtained from the first 
two moments. 

Table IV shows the mass and heat balance equations as 
applied to a tubular reactor (see Fig.l) having a laminar 
(Hagen-Poiseuille) flow profile given by 
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Ta^le_in 

Species Rate of F orm a tion and Heat Greneration Rate Terms 

= ki [C^][¥] - -2k2 [Pi]^o+21c2'[¥](?.^-[Pi]) 

-^3 iPiirc^] +>3* [P2] - k 3 [P 3 _][C 2 ] + k^* [P3] (b) 

H = ^1 - ki’[Pi] - ^2^0^ + ^ 2 ’ 

o 

+ [¥][C 2 ] - ^ 4 * [^2^ 

= k^ [C!2_][^] ~ k^ [P-|_] ■*" k^ “ k^ (Aq“ [P]_]) 

+ 2k5 [C 2 ] >0 - 2k5’ (A^ - [P^] - [P 2 ] ) 

+ 2k^ [¥][C 2 ] - 2 k 4 ' [P 2 ] ■ ■ (d) 

R^ = -k^ [C 2 ][¥] + k^' [P2] - k^ [02]?.^ + -[P^I-EPs]) 

(e) 

®¥ ' "’^1 + >^2^0 ■ - \'> 

-k [0 ][¥] + k ' [Pj] (f) 

4 ■ ^ ^ 

^ Q = Heat generation rate (Kcal/kg-hr) 

' ^ 

i=l ' 


( 6 ) 
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Table TII (C onti riued ) 


where , 


Rl 

= LCilLw] - k^ 

th] 


ih) 

R 2 

= ^2 ^0^ ■ ^ 2 * 

(Ai - 


(i) 


= 1^3 [0^] Tio - ^5* 

(Ao - 

[PJ) 

( 3 ) 


= ^4 [G2KW] - k^' 

tPj] 


(k) 

^5 

= ^5 [023X0 - 1 ^ 5 ' 

( Ao • 

- [Pi3 - [Pj]) 

( 1 ) 
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Grid pts 
forPO 
method 

^^Hypotheticol pt. 


(b) 

Fig.l 

( a) Schematic presentation of the tubular reactor with 
cooling jacket 

( b) Grid points for finite difference (methods I & II) 
and orthogonal collocation methods ( method IV) 
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(1-rVR^) (4) 

'Itie average residence titne# t, is given by 

t = (5) 

where z is the distance along the reactor. In writing 

these equations, the diffusion of all the components has 

been neglected and only thermal diffusion in the radial 

direction has been considered. While the diffusion of 

heavy polymer molecules may indeed be negligible, it may 

not be quite justified to neglect the diffusion of the 

lower molecular weight components, particularly, the monomer 

and water, This assumption is made for two reasons - 

firstly because precise data on the diffusion coefficients 

14 

are still not available , and secondly as the radial 
diffusion in the absence of vaporization is ejqjected to be 
a second order effect, and could be neglected at least for 
optimization studies. Tai et al. have also made similar 
approximations, it may be pointed out here that recently 
some work has been reported by Hamer and Ray^® on the 
detailed modeling of continuous tubular polymerization 
reactors incorporating radial diffusion of mass as well as 
property and velocity profile changes. Work along these 
lines for nylon-6 reactors could be carried out, but since 
our long range focus was on optimization studies, we did 
not pursue this approach. 



l3 


7 21 

The following correlations ' have been used for 
the various thermophysical properties 


thermal 

conductivity : k=0«21 Kcal/m-K-hr (6a) 


density 


specific 

heat 


zf 


+ 


(Kg/m^) = 1000^1.0065 + 0.0123 
(T(K)- 495) (0.00035 + 0.00007 


[Cl] 

[ClP} 


(6b) 


Kcal 


[ ^l] 


( = 0.6593 r + 


Kg-K 


[^1] 




^0.4861 + 0.000337 T(K)j 


[^] 1 X 

[ = l]o) 

(6C) 


The boundary condition at the wall of the reactor for 
the heat balance equation (eqn.h. Table IV) is obtained by 
considering a coolant at a constant temperature, Tj(*K)/ 
circulated in a jacket enclosing the reactor (Fig. la), the 
overall heat transfer coefficient being h. The symmetry 
condition at r=0 is given by eqn.i. Table IV. 

A radial temperature gradient will be established in 
the reactor at any cross-section along the length because 
of the presence of the velocity profile as well as the heat 
flux through the wall. As a consequence, the concentrations 
of all the species will also have a radial variation. Hence, 
the variables [ C^] , “Xq' ^l'[^2^ ' t Tare 

all functions of r and t in general, and can be represented 
as ] (r,t) , etc. Tb obtain all these profiles, the set 
of differential equations (eqns. a-g. Table IV) need to be 
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Table IV 


Mass and Heat Balance Equations in a Tubular Reactor with 

Laminar Flow Profile 


Mass Balance: 


2( 1 

2 /„2v 

- r /R ) 

dt 

(a) 

2( 1 

2/r,2v 

- r /R ) 

diPj 

dt ^1 

(b) 

2( 1 

2/„2> 

- r /R ) 

at 

(c) 

2( 1 

2/„2, 

- r /R ) 

dt ^1 

(d) 

2( 1 

2/„2. 

- r /R ) 

dt ^2 

(e) 

2( 1 

2/„2, 

- r /R ) 

d[ w] _ „ 

at ^ 

(f) 


Heat Balance: 

2(1 - rW)fC^ ||-= [r -f/CAQ) 

(g) 

(h) 


-k 

Sr 


= h (T 


r=R 


- 


St 

3 r 



(i) 


av 


t = z/v 
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integrated simultaneously with the boundary and symmetry 
conditions (eqns.h-i. Table IV). 

A close examination of the equations of Table IV 
alongwith the rate expressions of Table III reveals that 
the equations for [^2 1' ^1' t ^2^ ' ^ ’and T do not 

involve any other ^variable other than [ P 2 ] aJid [ P 3 ] . 

The equations for |■P 2 ] and [ P^] will similarly involve 
[P^] and [ P^] . Tto break this hierarchy of equations and 
to make thorn solvable/ a simplifying assumption/ called 
the closure approximation/ is employed. This is as follows; 

['’3] ' [^2] ' [hi 

1 22 

This was first used by Reimschuessel . Tai et al. have 
found that the final results are insensitive to this 
assumption. The equations of Table IV with eqn.7 now 
comprise a complete set of coupled differential equations 
which need to be integrated numerically. 

TWO techniques are employed in this work to convert 
the partial differential equation - PDE (eqn.g/ Table IV) 
into an ordinary differential equation - ODE. These are 
the finite difference (FD) and the orthogonal collocation (OC) 
methods. The basic principles of these methods are explained 
in Ref. 23 . 

Tt> use the FD method# the radial distance# 0 ^ r ^ R 
is first divided into N equal parts# of length Ar each# 
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using (N+l) grid points numbered from 1 to N+1 (Fig. lb). 
Ibus the jth radial grid point at is at r= ( j-1) ( A r ) . 

Ihe variables at r^ are renamed for the sake of convenience 
as follows: 

= [C^] (r^t) ; 5 . [P^] (r^t) ; 




= > 


0 


(r^^t) 


U4_.(t) . 


(r^.,t) 


u 


5, 



U .(t) 5 [W ] (r .,t) ; 
D/ 3 J 


U„ .(t) = T (r. t) ; j = 1,2, (8) 

' f 3 3* 

Ihe radial derivatives in eqn.g. Table IV, the 

boundary condition (eqn.h. Table IV) and the symmetry 

condition (eqn.i. Table IV) are converted into algebraic 

23 

expressions using the FD technique . This has beoi done 
in two different ways, in th e first method (method I), the 
Laplacian in the heat balance (eqn.g. Table IV) is broken 
up into two parts to give 


2(l-rVR^)/Cp|| 


r a r 


+ k 


^ + /( AQ) 
dr2 


( 9 ) 


23 

Ihe finite difference formulae for the first and second 
derivatives (mid-point formulae) are then applied to give 
following ODE corresponding to the jth point: 
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7 7 ^^7 • 

2(l-ryR') fjCCp) . 


= k [ 


U„ . . - U„ .1 
7/J+l 7,J-1 

2(A r) 


U_ . 1 - 2U„ . + U„ . T 

7,J+1 7/j 7^j-l 

r ) 2 


+ ; j=2, ,N (10) 

At the center, i.e., at j=l (r=0) eqn.lO presaits a problem, 

because it involves at an undefined point, j=0, as well 

as involves division by zero (r^=0). This is resolved by 
. $ 

using L'Ho^ital's rule as well as the symmetry condition in 
23 

its FD form : 


^ 1,2 “ ^ 7,0 
2(Ar) 


O; or, U 


7,0 


U, 


7,2 


( 11 ) 


Here, ^ is the temperature at a hypothetical point, j=0 
(Fig, lb). The resulting equation corresponding to j=l is. 


7 7 '^^'7 

2(l-ri/R ) 


4k 

(^r)^ 



^7,1> + 


( 12 ) 


Similarly eqn. 10 for j=N+l (r=R, at the wall) involves 
^7,N+2' temperature at another hypothetical point. This 
is resolved by using the FD form of the boundary condition 
(eqn.h. Table IV): 


^7,N+2 " ^7,N 


2(Ar)h , . 

k ^ 7,N+1 



(13) 
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Iliis is used in eqn»lO for j=N+l. The set of mass balance 
equations can be transformed to the FD form quite easily. 
Thus, 7(N+1) ordinary differential equations have been 
generated, which are given in Table V. The number of 

variables to be solved for are also 7(N+1). These DDEs 

23 

are solved using a GEAR'S package , which has a built-in 
step control algorithm and is particularly useful for stiff 
systems. The NAG Library Routine D02 ebf was used for this 
purpose in a DEC 1090 system. 


The second approach of the finite difference technique 
(method II) is very similar to method I except at the 
boundary points, j=l and N+1. To apply the FD technique to 
the boundary and symmetry condition, one-sided first 
derivative formulae (ref, 23, p.68), which have the same 
degree of accuracy as the mid-point formulae used in method I, 
have been used. The symmetry condition at r=0 is transformed 


into: 



4 

3 




(14) 


and the boundary condition at r=R gives. 


u, 


7, N+1 


~ 2 (A r) ^^7,N-1 “ '^^7,N^ 


h + 


3k 


(15) 


2(Ar) 


Use of the FD form of the heat balance equation at the 
center and boundary points is avoided by using eqn. 14 and 15. 
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Table V 


Set of Differential Equations After Applying Method I Of 

Finite Difference Technique 


2(1 - r^/R^) 


2(1 - r^/R^) 


2(1 - r^/R^) 


2(1 - rJ/R2)J>i,(Cp)^ +/^(AQ)^ 


du. 

i 

dt 

du_ 

^ i 

dt 


dt 


6J 


*'\V 

j = 1, — 

.... N+l 

( a) 


J !/•••• 

. . . .N+1 

(b) 

: 





j 

M+l 

(c) 


dU, 


7,1 _ 4k 


dU. 




(d) 


= kl 


r 1 


^7, j+1 ~ ^ ^7,1+1 ~ ^^7,j ^7,3-1 


•j 


2( A r) 


] 


( ^r)‘ 


+ ^ j = 2,. 


,N+1 


with 


^7,N+2 ^7,N 


li^flJ3(u 


“7,N+1 ” 


( e) 


dependence of variables on t is not explicitly mentioned 
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Ihe resultant set of odes consists of eqn. a-c. Table V 
for j=l,.../N+l and eqn.e. Table V for j=2,,...,N. Ibe 
total number of Odes is 7(N+l)-2 and the number of 
variables is also 7(N+1)“2 , The other two variables/ 

i.e./ the temperatures at the center and at the external 
boundary point, can be obbined from eqn. 14 and 15. Ihe 
same GEAR'S package is used for integrating the set of 
equations as well. 

A compact computer package for integrating PDEs (NAG 
Library Routine D03PBF) was also used (method III). This 
package consists of a combination of the finite difference 
technique and the GEAR'S package to solve a set of PDEs with 
appropriate boundary conditions, as discussed later, 
method I and III are essentially the same. However, method III 
has an additional flexibility in the sense that it can be 
used for any type of unequal grid spacing, upon changing an 
input integer, imesh, Ihe value used for our case was 4, 
which results in having more grid points near the external 
boundary (where the gradients are expected to be steeper) 
according to formula: 

r^ = R Sin ( ) (16) 

Ihis method was used to generate various simulation results. 

Another technique tried was the method of orthogonal 
collocation. The ability of this techniques to give results 



of greater accuracy with less computational effort under 

23 

certain conditions, is well established . However, no 
work has been reported till now employing this efficient 
method for the simulation of nylon-6 tubular reactors 
with radial gradients. In this method (method IV), the 
variables are approximated as a series of orthogonal 
Legendre polynomials. Hiese polynomials are symmetric 
in nature with respect to the independent variable which 
is X = r/R. Ihe polynomials thenselves are obtained from 
the orthogonality conditions e.g. , 

1 

W(x^) (x^) (x^) dx=0, 

° k ^ m-1 (17) 

2 

where, W(x ) is a wei^tage function which has been taken 
2 

as ( 1-x ) in the present case. 'a' is the geometry factor 

and has the value of 2 for the cylindrical geometry. Once 
2 

W(x ) and ‘a* are fixed, all these polynomials can be 

obtained. The polynomial with the highest degree to be 
S ' 2 

considered ijt p (x ) where N is the number of interior 

N 

collocation points (excluding the external boundary). The 

V * 2 

roots of Pjg(x ) are taken as the interior collocation points. 
Now the matrices A. . and B, . may be formed^^ from the 
collocation points so that, for any dependent variable y. 


dx 


at ith point 


N+1 

i: 

j=i 




(18) 
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and 



at ith point 


N +1 


i: 

■ 5=1 


B. . 
1/ J 


(19) 


using eqn.18 in the boundary condition (eqn.h. Table IV) 
for our case, after nondimension ali zing the independent 


variable r in terms of x, we get< 

n, N 


U 


7,N+1 


R ^+ 1,3 


+ hT, 


( 20 ) 


^ R N + 1 


ThuS/ there is no need to apply the heat balance at the 
boundary. Ihe symmetry condition need not be considered 
here as, this is automatically taken into account by the use 
of the symmetric polynomials. The heat balance (eqn.g/ Table IV) 
after being written at the jth collocation point gives. 


o dU , 

2(1-Xj) / dt'^ 


N+1 


= ^ XI ®j,i ^7,i 

j=l,....N (21) 


i=l 


After substituting eqn. 20 in eqn. 21 the final form of the 
heat balance isi 



9 ^ : 

zn-.]) 


N 


( — ) 


B. . U . 4 
3# i 7,1 


i=l 


N 




^ r2^ ®j,N4l 


(h 4 


^41, 


) 

+ /.(AQ) 

J J' 


N4l 


j=l,....,N (22) 


Ihe mass balance equations (eqn, a-f. Table IV) can 
again be applied at all the (N 4 I) collocation points, Ihe 
resultant set of [ 7(N4l)-l]0DEs are listed in Table vi. Frcwi 
these, the [ 7(N4l)-l ] variables (other than N+l' 
temperature at the boiandary ) can be obtained by the use of 
GEAR’S method with appropriate initial conditions. The 
temperature at the boundary (u^ N 4 l^ obtained by using 
eqn. 20. 


In all the numerical techniques, the variables obtained 
at the radial grid points at any t aieused to -get the average 
values. The averages can be computed using the 'cup mixing 
rule' to give, for any concentration and moment. 


R 




j' r(l-rW)Uj^/(r) dr 


av 


J r ( 1-r VR^) / (r ) dr 
for i“l, . • . , 6 


o 


( 23 ) 
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Table VI 


Set of ODES 


2(1 - 


2(1 - 


dU.. .. 


dt^ = .3 = 1 N+1 






= 1 , 


2( 1 - ^ j ) 


. j = 1. 


2(1 -X.") /■j(Cp)j ^ X 


*p2 ' ®J,N+1 


rI Vl,i “7,i * >' Tj) 

1 = 1 


( h + — a 3 

R ^-N+l, NH*1^ 


+ f j(^Q)j ^ 


3 = 1 . 
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and for temperature. 




J r(l-r-?'E^U7/(r)Cp(r)ar 


R 

j" r(l-r}R^f(r) Cp(r) dr 


( 24 ) 


In the first three methods the integrals are computed 
by anploying a third order finite difference formula given by 
Gill and Miller^^. Ihe NAG Library Routine DOlGAF is used for 
this purpose. In the orthogonal collocation method, the 
integrals are evaluated using the weightage vector generated 
from orthogonal polynomials . Ihe final form is as followss 


for concentration and moments. 


(up 


av 


J — , 

N+l 


i~l, ... • , 6 < 


( 25 ) 


and, for temperature, 

N+l 


(U^) 


av 


N+l 

X W, (l-X^.) /, (C^) 


3=1 


j 


r ■' j ' p' j 


( 26 ) 


In laminar flow, the residence time is infinity at the 

-4 

wall. To avoid numerical problems, a small slip (v^(r=R)/’\^^=lG ) 

was assumed at the wall. Ihe nmerical results were found to 
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be relatively unaffected by changing the dimensionless 
velocity around 10 Ihe inlet condition (t=0) at the 
wall (r=R) was also modified. The equilibrium conversion 
values of isothermal runs at 250®C were fed in as values 
at r=R/ t=0. This reduced the computational time significantly 
and avoided convergence problans. Any change in these input 
values Was found not to affect the results, but led to longer 
computational times. 



CHAPTER 3 


RESULTS AND DISCUSSION 

A general computer program was made. Results for 
isothermal conditions could be obtained by putting aq=o in the 
computer program. Similarly, adiabatic operation could be 
simulated by putting the heat transfer coefficient, h, as 0. 

The program could be run using either plug flow (v = Const ant= 
or parabolic flow (eqn,4) profiles, A typical run for 
t = 20 hrs and 10 FD points took .3 min 45 sec and for 10 OC points 
it took 2 min "36 sec on a DEC 1090 system. 

The computer programs for all the four methods described 
previously were checked with results available in the liter at ure, ^ 
The programs for methods I— IV were run for the simple case of an 
isothermal plug flow reactor with the input conditions as* 

[C^] ^ = 8,8 moles/kg 
[W] ^ = 0,16 moles/kg 

[^l] o =(>> 0)0 = <^^o = [^2] o-«-- ^27) 

and for tenperatures of 230‘’C, 240'C, 250“c, 260®C and 270*0. 

The results were found to be in complete agreement with those 
25 

of Ray and Gupta, thus confirming the correctness of the computer 
programs( except the heat balance equations). The programs were 
also run for an adiabatic plug flow reactor with the feed coming 
from the top of a VK column^^. The results were found to match 
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those of Gupta and Gandhi^*^. This further conf firmed the 
correctness of our computer programs { including the heat balance 
equations) . 

Having developed confidence in the computer packages 

prepared, detailed temperature and concentration profiles for 

several conditions were generated. The feed conditions were 

the same as in eqn.27. Table VII gives the various flow and 

heat transfer conditions used. Figs. 2 to 17 present profiles 

for monomer conversion, number average chain length(u.^), 

tenperature and cyclic dimer concentration at various values 

of t (S 2 /v^^) , For heat transfer coefficient, h, a typical 
11 ‘ 2 . 

value of 5,0 kcal/m -K-hr was used. The existence of distinct 
concentration and temperature profiles are quite evident from 
these figures. 

Figs. 2 to 5 show the profiles at t = 5hr for various flow 
and heat transfer conditions (Table VII). The laminar flow case 
( ^curves & 6) is the most Interesting, Since the velocity 
near the center is higher, the residence time of a fluid element 
near the center ( at the same t) is smaller when compared to that 
of a fluid element near the wall. Thus the reaction is close 
to equilibrium at the wall. A maximum in the temperature slightly 
away from the wall is observed ( curves 5 & 6, Fig-2) at t=5hr. 

This represents a combination of higher heat gene artion ( due to 
higher near equilibrium conversion) near the wall as wSll as heat 
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Fig.l3. 

Radial profile for at t=15 hr for various flow 

and heat transfer conditions (given in Table VII). Feed 
condition is as per eqn. 27. 




Radial profile for temperature at t=20 hr for 
various flow and heat transfer conditions' ( given 
in Table VII). Feed condition is as per eqn. 27 








Radial profile for at t=20 hr for various flow 

and heat transfer conditions < given in Table VII). 
Feed condition is as per eqn. 27. 
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transfer effects. In the laminar adiabatic case ( curve 5, 

Fig.2), a temperature maximum first appears ( at t near zero) at 

the wall because of almost instantaneous attainment of equilibrium 

conversions there. The higher temperature leads to thermal 

diffusion away from the wall as t increases (since h=0), contri- 

bgting to a decrease in the temperature at the wall. In addition# 

conversion 

as t increases layers adjacent to the wall reach near-equilibrium/ 
(curves 5 & 6, Figs. 3 to 5)# generating more heat there. This 
explains the temperature maximum in curve 5, Fig.2. The peak 
shifts further inwards as t increases (fig. 6). Because of 
cooling from the jacket# the temperatures at the wall are lower 
for the laminar non-adiabatic case (curve 6# Fig.2). The slightly 
higher near-equilibrium conversions for curve 6# Fig.2 (compared 
to cuirve 5), is because of the lower temperatures present and 
the overall exothermic nature of the reactions. Since the 
equilibrium values of are Bore sensitive to temperature# the 
difference between curves 5 and 6 (Fig, 4) are greater. In plug 
flow# the residence times of fluid element at the wall are not 
infinite# and so# near-equilibrium conditions do not prevail at 
the wall for this ideal reactor. In the adiabatic case# the 
temperature ( and other) profiles are flat as expected# since 
there is no driving force for creating radial difference®. For 
the plug flow non-adiabatic case (curve 3# Fig.2) a temperature 
maximum near the wall develops again. Near t=0# a maximum 
appears at the wall (Tj^T^). This beds to higter conversions 
at the wall# which simultaneously leads to higher heat genera- 
tion by reaction at that point. Very soon the temperature 
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at the wall increases beyond T . Heat transfer takes place 

U 

both to the jacket as well as towards the center. The conversion 

is the highest at the wall for case 3 (Fig.j) since it represents 

a cumulative effect from t=0 to 5. Figs, 6 to 17 show conditions 

at t=10, 15 and 20 hr. The trends at t=5 hr continue, but some 

flattening of temperature profiles in observed, because of thermal 

diffusion effects away from the peak, both towards the wall as 

well as towards the center. The conversion, |Jt and cyclic dimer 

n 

profiles show higher levels near the center as t increases, as 
expected. Figs. 15 and 17, at t=20hr, show that even though the 
monomer conversion has almost attained equilibrium values throughout 
the cross-section, lower cyclic dimer concentrations near the 
center are observed. This is because equilibrium for the cyclic 
dimer reactions has not yet been attained. 

Figs. 18 and 19 show the profiles of average temperature, 

T, and average chain length, u^, along the reactor length for 
various flow and heat transfer conditions. It is interesting to 
note that the distinct "S* - shape for the plug flow case (curves 
2 & 3, Fig, 18) is relatively less prominent in the laminar flow 
case (curves 4 & 5, Fig. 18). Average temperatures for the 
non-adiabatic ease (curves 3 & 6, Fig, 18) are in general lower 
than that for the adiabatic case (curves 2 & 3, Fig. 18) because 
of the cooling effects of the jacket fluid, except for the initial 
part of the plug flow case (curves 2 & J) where the jacket fluid 
has a heating effect as discussed earlier. For 'jl^, again, the 
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conditions (given in Table VII) with feed condition as 
eqn. 27. 
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'S' - shape in case of plug flow ( c'.irves 1,2 & 3, Fig. 19) is 
less prominent in laminar flow (curves 4,5 & 6 Fig. 19). In 
general, isothermal operation (curves 1 & 4, Fig. 19) gives 
lower in the initial part of the reactor becaBse of comp- 
aratively lower temperatures in this case. But as t increases 
and equilibrium conversions are attained, for the isothermal 
case registers higher values. 

' Table VIII shows the effects of the various parameters, 
e.g., radius of the reactor (R), coolant temperature (T_), feed 
temperature (T ) and inlet water concentration ([w^ ), on the 
average values of temperature, conversion, |i ^ and cyclic dimer 
concentration. Decreasing the value of R to 0.4 (Run-2) from 
the reference conditions (Run-1) helps increase radial heat 
transfer. Thus, the initial heating effects of the jacket 
fluid Increases T at 5 hrs <0ver that for the reference run. As 
t increases, the temperature is comparatively lower due to 
better cooling. However, average conversion and ]ij^ are always 
lower. Increasing the value of R (Run-3, Table VIII) has the 
opposite effects, A lower coolant temperature (Run-4, Table 
VIII) causes lower temperature in the reactox. This reduces 
the average conversion and initially. But, as equilibrium 
is attained at higher values of t, the average conversion 
and u^, are higher as an effect of lower temperatures. A lower 
feed temperature ( Run-5, Table VIII) causes lower T. As 
before, average conversion and are initially lower, but 
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I !. T.. KANPUP. 

Table VIII A.iOM?,l 

Effect of Various Parameters on Average Values 


Run Description 

t(hr) 

T , 
(“C) 

Avg.Conv. , 

% 

[C^] 

1. Reference Condition: 

5 

264.09 

28.24 

67.08 

0.00516 

Laminar flow. 

Non- adiabatic (h=5) 

10 

287.20 

73.66 

123.83 

0.01788 

R = 0.6m 

T = 260 *C 

15 

294.09 

88.50 

149.61 

0.03003 

T^= 250 ®C 

Feed condition as 

20 

293.83 

89.65 

151.36 

0.03702 

per eqn.27 






2, R^ = 0.4m 

5 

264.27 

29.09 

68.79 

0.00497 


10 

287.06 

74.89 

126.46 

0.01782 


15 

292.90 

88.80 

150.72 

0.03002 


20 

292.01 

89.78 

152.52 

0.03684 

3. R^ = 0,8m 

5 

263.96 

27.66 

65.69 

0.00524 


10 

287.39 

73.19 

122.80 

0.01786 


15 

294.71 

88.38 

149 .09 

0.03000 


20 

294.73 

89.59 

150.80 

0.03706 

4. = 250 ®C 

U 

5 

263.49 

28.03 

66.71 

0.00486 


10 

286.44 

73.67 

124. 19 

0.01751 


15 

29 3.12 

88.55 

150.22 

0.02969 


20 

292.64 

89.73 

152.13 

0.03670 

5. T® = 240*C 

o 

5 

247.96 

15.38 

53.06 

0.00253 


10 

264.57 

47.73 

89.55 

0.00940 


15 

278.46 

74.90 

128.58 

0.01757 


20 

284.75 

88.09 

153.77 

0.02685 
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Table VIII (Contd.) 


6- ["lo 

= 0.13 moles/kg 

5 

260.25 

20,76 

65,07 

0.00352 



10 

280.21 

60.04 

117.87 

0.01301 



15 

292.05 

84.16 

161.87 

0.02413 



20 

29 3.84 

89.17 

170.77 

0.03260 

7* 

=0.18 moles/kg 

5 

266.66 

33.24 

67.38 

0.00638 



10 

290.52 

80.16 

125.50 

0.02108 



15 

294,38 

89.30 

140.50 

0.03305 



20 

293.72 

89.71 

141.25 

0.03892 


^All otber value as in reference condition. 
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are higher near equilibrium. Lower feed water concentration 
(Run. 6, Table VIII) causes lower conversion and]!^ initially, as 
the ring opening step (step.l. Table I) is slower. But as t 
Increases, and near-equilibrium conditions are attained, 
conversion improves and reaches the values for the reference 
state, since equilibrium conversion of polycondensation reaction 
{step, 2, Table I) is higher for lower [ w] , But the effect 
on is much more pronounced. This is because a shift of the 
equilibrium for the polycondensation step to the right hand 
side increases the chain length considerably. The effects of 
higher initial water ccMncentration (Run. 7, Table VIII), are opposite 
in comparison to run. 6, Table VIII. 

The effect of some 'computational variables' is now 
studied. The values of the 'slip', i.e., the non-dimensional 
velocity ( v^ (r=R)/v^^) assumed at the w«ll(to avoid overflow and 

other computational problems) was changed aroung the reference 

-4 

value of 10 and the effects are shown in Table IX. It is 
observed that the results are relatively unaffected by changing 
the slip, but the CPU time requirement increases for lower values 
of slip. Similarly, instead of using the actual feed conditions 
at t=0, r=R, values of t^lJo' close to 

equilibrium values at the wall temperature were used, for the 
laminar flow case, to makes the equations less stiff. These 
values were changed to see their effect on the numerical results 
further down the reactor. The results were found to be quite 



Table IX 


Effect of Slip on Simulation Runs 


Run^ 

Slip used 

Wall Temp, 
at 5 hr (®C) 

CPU 

Time (for t= 5 hr) 

1 

lO^^Cref ) 

287.51 

2 min 25 sec 

2 

10“^ 

287.51 

3 min 5 sec 

3 


287.51 

2 min 6 sec 


^All runs are for laminar flow adiabatic rector with feed 


condition as per eqn. 27 
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insensitive, but feeding in the near-equilibrium values reduced 
the computational time. 

Fig. 20 shows the effect of increasing the number of FD 
points on the results at t=5 hrs. The runs were made using 
method III, with the spacing of grid points denser near the wall 
(eqn.l6). The figure' shows very clearly that the results, 
particularly near the wall, converge as the number of points are 
increased and after about 8 points the results do not show any 
major discrepancy. 

Methods I and III (for identical placement of FD points)wsre 
found to give the same results for various reactor conditions. 
This confirms that these two computer programs, one in the NAG 
Library and one made by us, are essentially the same. Both of 
them use the FD technique with GEAR'S method. Method III, 
however, is more flexible since it permits the use of unequal 
spacing of FD grid points. 

Some discrepancy near the wall was found between methods 
I and II. Table X shows the wall temperatures (which is the most 
sensitive variable )fur the laminar flow adiabatic reactor as 
predicted by methods I and II, at t=5 hrs. With 9FD grid points, 
the results show some difference$at the wall, though those at the 
internal positions agree well. The reason for this discrepancy 
lies in the fact that in method II, the wall temperature is 
obtained using the boundary condition alone (eqn.l5). Therefore 



t = 5hr 
— lOpt PD 
o 4pt FD 
* 6ptFD 
o 8ptFD 


Adiabatic reactor 
Laminar flow 


X X 


Fig. 20. 

Temperature profile at t=5 hr in a laminar flow adiabatic 
reactor (feed conditions given in eqn. 27) for different 
number of grid points using finite difference (method III) 
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Table X 


Comparison of Methods I and II 


( Laminar Flow Adiabatic Reactor) 


No. Of 

Temp, at 

O 

the Wall ( ) / C/ at t=5hr 

r U pu S • # "" 

( N+l) 

Met hod -I 

Method-II 

9 

290.04 

294.85 

12 

285.72 

288.68 

15 

287.62 

288.04 

20 

287.50 

287.50 
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the ODE from the heat balance equation at the wall (r=R)/ and 
thus the heat of reaction information near that location, is 
not used. In method I, on the other hand, both the boundary 
condition and the ODE from the heat balance are used at r=R. 
Therefore, when the gradients are steep near the wall ( for lower 
values of t), method I gives better predictions. Thus method I 
is superior to. method II, particularly with regard to the 
results near the wall. However, as the numbers of FD points 
are increased, the results for both the method converge as 
seen in Table X. Similar results are obtained for the non- 
adiabatic case. However, for the non-adiabatic case of the 
flow reactor, where the temperature gradients near the wall 
are not so steep, the results agree well for the two methods 
( 1 & II) . 

The orthogonal collocation technique gives results which 
match those from method III for all cases of plug flow reactor 
( run. 1-3, Table VTI) as also for the isothermal laminar flow 
reactor (run. 4, Table VII ). However, the predictions tio not 
match well for the adiabatic (see Fig. 21) and non-adiabatic 
cases of the laminar flow reactor (Runs. 2 & 3, Table VII). 

The mismatch is particularly large for lower values of t, when 
a steep gradient of temperature exists near the wall. The 
ability of our computer program for the OC method to predict 
good results for cases where the temperature gradients are not 
too steep near the wall, indicates the correctness of the program. 



59 



Fig. 21. 

Temperature profile at t=5 hr in a laminar flow adiabatic 
reactor (feed condition given by egn* 27) for different 
number of grid points using orthogonal collocation (method IV) 
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It is felt that the incorrect results obtained for the 
non-isothermal cases of the laminar flow reactor at lower values 
of t could be attributed to the fact that like method II, the 
OC method also does not incorporate the ODE from the heat 
balance at the wall. Instead, it computes the wall temperature 
(Tn^i) using all the other values of temperature, 

( eqn.20). The OC method does not compare too well with method II 
either. This is due to the fact that in calculating the wall 
temperature method II uses information only for the two 

adjacent points, T„ and T„ - (eqn.15), whereas the OC roe thud 
requires information ef all other points (eqn,20), giving 
different weightages to them. Thus the predict^ion of the wall 
temperature in the OC method is likely to be sensitive to the 
position of the maximum { as it shifts) in the temperature profile. 
As t increases, the temperature maximum flattens out and the 
OC method is indeed found to give better results. The results 
for the OC method, where steep temperature gradients exist, are 
likely to improve on the introduction of the finite element 
technique. 



CONCLUSIONS 


A detailed simulation of nylon-6 polymerization in 
continuous flow tubular reactors was carried out using four 
methods. Two types of finite difference techniques as well 
as the method of orthogonal collocation were used to convert 
PDEs into ODES. GEAR'S method was used to solve the 
resultant set of odes. 

The simulation results showed the existence of 
significant radial gradients for temperature and concen- 
trations for the laminar flow cases, particularly for lower 
values of t (average residence time). A temperature peak 
forms at lower t near the wall, which flattens out and 
simultaneously shifts towards the center further downstream. 
Ihe monomer conversion and number average chain length 
profiles result from the temperature profiles. Higher 
temperature normally increases conversion and number 
average chain length, but it has a reverse effect where 
near- equilibrium conditions exist, a characteristic *S'- 
shape was obseirved in case of plug flow in the plot of 
average temperature and versus t. The plot was much 

flatter in the laminar flow case. The cases of adiabatic 
and non-adiabatic reactors both with plug and laminar flow 
were found not to give much different results for the average 

value of , which is of greater importance as far as the 
n 

final product properties are concerned. 
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Methods I and III were found to be essentially 
similar with the latter having more flexibility as it 
could use any type of spacing of grid points. Method II 
was found to give poorer results near the wall particularly 
when the temperature gradient near the wall was steep. 

Both the methods i and II agreed wall as the number of 
grid points was increased. Method IV gave good resxilts 
for all the cases except the adiabatic and. non-adiabatic 
cases in laminar flow situation near the wall where the 
temperature gradients were particularly steep. Use of 
the finite element method is suggested to give improved 
residts for such cases. 
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FINITE DIFFERENCE! METHOD I 

*^*:^*:$*******t**^^*t*=¥^*$***t******^********^* 
RS!1§K§^9n„S%^?l^i^P^®^'“^150,200),RX(25 3 

common /OUT/H»IT,a»N 

COMMON ICOL,COL£U10,RX,UO#TW,R,TK,U70,ICQND,TVEL,DR 

external FCft,OUtPUT,REDERV 

open C UN I Ts; 15 , DEVICES 'DSK % FILES ' FD . IN ' ) 

OPEN CUNITs 39,DEVIC|=»DSK*, FILE* ^FD. OUT 
OPEN CUNITse, DEVICES 'pSK', file* 'ERR2 ,001*) 

READa5,») IC0ND,IVEL,WPTS, 

1 UO,TK,TW, 

2 R aENo 

3 Ul6,U20aj30,U40,U50,U60,U70 
IC0L=NPTS»1 

IPCIC0ND.EQ.D3 WRITE{39,17) 

IFCICOND.EO.l) WRITEC39,18) 

IF(ICONO,E0.2) «RlTE(3<}/19) 

FORMATC15X,50C '♦'),///, 30X, 'ISOTHERMAL RFACTOH ',///, 15X ,50 ('*' ) 

I'OR^AT(15X#50('»'),///,30X,' ADIABATIC REACTOR ',///, 15X, 50 ( '*' V, 

loRMAici5X,50C'»'),///,22X. 'TUBULAR REACTOR NJTH HEAT 
1 transfer '»///,15X.50 ('**),//) 

WRITE(39,203 R,NPTS,ai0,U40,U20rUSOj,U70,U30,U60 
FORMAT(5X,'RA0lUSs',F6.2,5X,'NO OF F.p. P0TNTS=',I3,/, 

1 5X# 'MONOMERS', r8,4,5X^'lst MOMENTS' fF8 .4 ,5X, / , 

2 8X,'P(ns',F8,4.!0X,'DIMER='#F8.4.10X,'TEMP= ',F10.4./, 

3 2X»'0th MOMENTS*, F8. 4, lOX, 'WATERS*, F8, 4,/, 75('i')) 

TN=TW-273. 

IFCICOflD.EQ.2) WRITE(39,213 UO,TN 

FORMAT(/,$X, 'HEAT TRANsF COEFFs ' ,F10.5,2X, 'KCAL/f SO M)(HR) 

1 (DEG K) 'COOLANT TEMP= ' ,F10,5 , // ,75 ( '♦ h 3 

TOLsI,E-06 

IRELA§=2 

MPED=0 

XsO.O 

HsO.5 

IT=0 

DRsl.O/lCOL 

DO 11 Isl,ClC0L+l3 

C0Lh)=FL6AfCI-l)*DR 

DO 33 Jsl.ICOLtl 

RXCj)s2,*h.-(COLCiJn*»2.3 

IFCRXCJ3.E0.O,O) RX(J)sl,E-04 

IFhVEL.EQ.6) RXCO)s|,0 

continue 

initial condition 

DO 22 Isl,(lCQL+l) 

Y(I3=U10 
YCICOL+1+I)=U20 
YaC0L*2 + 2 + I)sU30 

y(IC0L»3+3+I)=U40 
y(IC0L»4+4+I)sU50 
YaC0L^5+5 + l5=U60 
Y(ICOL*6-I’6+I)5U70 

rfiiffiii# 

YCIC0L*4+4)=8, 05941 
yaC0L*5+5>=0. 034956 
Y(ICOL»6+6)=0, 116415 
IFAILsO 
fIsClC0L+l>^7 
Ili=ia+N 

CALL D02EBF(X,XEND,N,Y,I0L,IRELAB,FCN,MPEB, 



6B 


26 


24 

C 

C 

C 


44 

C 


16 


11 


,EC(5},C0L(25},HHUf2f!j 


m pederv.output.w.iw.ifail) 

IFdFAlL.ME.O) GO TO 24 
f|pE »,IP’AIL 

1 1 i » f » ♦ ♦ ♦ f f t f 4- ^ )j[ jf; ^ ^ 3|5 ^ 5,; ;|S ^ g ^ ^ 3(. ;|t I(t 3|: ♦ !jt f f 4! ♦ * ^ # 

SUBROUTINE FCNCTdrF) 

DIMENSION YC15d);F(150),AO(5),E0(5),FK(5,25) 

1 ®!Ji^54i?5APj[5<25)jDS(5),DH|5),ACf'- ' 

2 ,CP(25),RXC2S5fTP{253,FlC25),F2C2..^ 

COMMON lCOt,COli.UlO,RX,u6,TW,R,TK,U70,ICOND,TVEL,DR 
DO 44 Jsl-lCOL+1 

TPfJ)=y(lCOU*6+6+J) 

TPC|COL+2)=TPCICOL)-2.0»DR»UO*R/TK*CTP(ICOL+13-TW) 

TYPE »,TPCICOL),TP(ICOL+I),fP(ICOL+2) 

DO 16 a=l,lCOL+i 

?HOC JJfiOOO. 0/Cl. 0065+0. 0i23*YCJ)+CTPCJ)-495.)* 

3 (0. 00035+0, 00607*YCa))) 

CPCJ)=0.6593*YCJ)/U10+{1.-YCJ)/U10)*(486.1+0.337*TP(J) 

4 )/lO0O.O 
CONTINUE^ 

A0(l)=5.9874E+05 

A0(2)=1.8942E+10 

AO(33*2,8558E+r?9 

AOt4>=8.5778E+ll 

AOC5>=2,5701S+O8 

EOCl)=1.988E+04 

EOC2/=2.3271E+04 

E0(3)*2,2845E+04 

iOt4>=4,2E+04 

E0(5352,i3E+04 

DH(1)=1.918E+03 

DHC2)=-5.9458E+03 

DH(3)=-4,0438E+03 

DHC4)=-t.6E+03 

DHC5)=-3.169lE+03 

DSCl)=;-7.8846 

DS(2)=9,4374E-Q1 

DSC3>s-§.9487 

OS(4)=-14,52 

DS(5)55.8565E-01 



• m 'm ♦ #■# '• # • 



,i)=l,§515E+07 
!)3t.2lHE+iO 
i)5l.i377E+|Q 
.3307E+I2 
,0UE+09 
_ .8806E+04 
i=2.067E+04 
S2.0107E+04 
?3.74E+04 
(a2.04E+04 
,987 
1=1 f 5, 

v'u Ai 0— l.ICOli+l 
FK(I.J)sAO(I)»EXPC-EO(I)/(UR*TP(J)))+ACCI)»EXP 
I (-ECCI3/CUR*TPCJ)))tY(ICOL»2+2+J) 
EK(I,a)=EXPCDSCI3/UR-DU(I)/CUR»TPCJ))) 
BKCl,0)=FK(I,J)/EK(I,J) 

CONTINUE 

DO 12 J=1,IC0L+1 

P2=YCIC0L+1+J? 

P3=Y(IC0I*+l+0) 

APRsCYCICDD*2+2+J)-2.*YtICOD+l+J)) 

IF ( APR . UT . 0.0) APR=0 , 0 
IFCAPR.DT.O.O) P250.O 
IFCAPR.DT.O.O) P3=0.0 
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C 

12 


13 

C 


15 

C 


17 

C 


f?!a6f;;i;j5: 

1 •»2,*FKC2. 

2 +2-*BK(; 

3 "FKprJ 

4 BKp, J 3 »P2-rKC5,JUyClCOIi+l+a)»YCICOL*4 + 4+J)+BK(5,J') ♦Pi 
FCICOL+l+J)aFCICOL+l+J)/PX(J) 



F(IC§Lf|+§+jJafK?l,a5JfU5*y(ICOL*5 + 5+J)«BKli,J)^y?iaol.+ l+Jj 
1 -FKC2#J5fKlCOL*2+2+a3»Y€lcaL'^2+2+J) + 

BKC2«J3?1fClCOIi»S + 5+J)*CYClCOL*3 + 3+d)-Y(TCnL,*2+2+J)} 


+FKC4#53JYflCOli»4+4+a)*y(ICOB^5+5+J)- 
BKC4,J)*P2 

COL*2+2+J)-FClCOL»2+24d)/RX(J) 


3 

4 

FCICOli* 

^cic6C^I+§+53sfKci#a3HcJ3lY?icoL♦§+§+j5-S^^H^n♦y(!auL^-l+.TJ 

1 +FKC3,J)#yCJ)»y(IC0L»2+2+J)- 

2 BK(3,J>*Cy(lCOB*2+2+d)-yCTCOL+l+J) J+2.»FK(5,J)* 

3 YCICOL4t4+4+J3*YClCOL»2 + 2+J)-2,*BK(5,J)^ 

4 APR+2.*FK(4,d3*YCIC0L»5+5+J)*Y(TCnLi*4 + 4+Jl-2.*BKC4,J)*P2 
FCICOL*3+3+J)=FClCOL*3+3+J)/RX(J) 


#(ICOi:i»4+4+J3=-FK(l#J)*Y?iaOL»4 + 4+J)»Y(ICOI.»5+5+J)+BKU,J)»p? 

1 -FKCSf J)*YCICOL*4+4+J3»Y{ICOL*2+2+J)+ 

2 BKCS-vDYAPR 

FC1COIj*44-4+J)=F(ICOL*4+4+J)/RX(J) 

F!icoBJ§+§+j3=-FK(n33ncd)^Hic6L^5+s+j34SKf ilaHYucoirfi+d ) 


^ FKC2,J)*YClC0B»2 + 24-J)fyaC0Ii=«=2+2+J)- 

2 BK(2,J)YYCICQ1.45 + 5+J)*CY(ICOL*3 + 3+J)-Y(ICOL»2 + 2 + J) ) 

3 -FKC4,J)*YClCOL»4+4+J)*y(ICnL*5+5+J)+RK(4,J) 

2 ♦P2 

FCICOLf5->-5+J)=:F(lCOL*5+5+vJ3/RX(J) 

CONfiBUE”'**’**’****’'”**”**’* 

DO 13 JJsi.ICOL+l 
P2=YCICOL+l+JvI) 

P3=Y{IC0L+1+Ja3 

APR=CYClCOL*2+2+JJ)-2.»Y{ICOL-l't+Jj)) 

IFUPR.Lj^g.L) 

IFCAPR.LT.O.OJ P2=0,d 
IF(ftPR.liT.0,05 pSsO.O 

ri(JJ)s-DH(i)»CFR(l,JJ)*YCjJ)»y{ICOL*5+5+JJ)-aK(t 
7 yClCOL+ltJJ3) _ 

-DH(2)»CFKC2,JJ>»Y(TCDL^2+2+JvI)»Y(XCOL*2t2f JJ)“Btv(2,Jd)» 
YCICQli*S+S+aj)»CY(lCOI.^3+3+v.W)-y(ICOL^2 + 2 + JJ)) j 
-DHC33»CFKC3,JjpY(jJ)*y(TCnL*2 + 2+J*J)-flK(3,JJ)4= 
YClCOL*2+2+Jd)+.BK(3,JJ)*Y(ICQB+l+JJ)) 

-DaC4T*CFKf4,JJ)n(ICCL*4+4+vJd)*Y(ICOL»B+5+JJ)-0.\(4,id>*P2 • 
DHC5)*CFK(5,JJ54Y(IC0L44+4+dJ)*Y(IC0B<=2 + 2+JJ)*B6C5,j.t)fAr'.O 


3 

4 

5 

6 


CONTINIJE 


DO 15 L=2,ICOL+l 

F2(L) = (TK/R*42,(n*((7P(L+l)-TP{L-l))/2.0/l)R 

1 /COL (li ) + ( TP ( L+ 1 ) -2 . 0»TP ( L3 +TP C L-1 ) ) / C nR*^2 . 0 ) ) 
COnilNUE 

TYPE ♦.Fi(ICOL+n,F2(lCOt+l) 
F2(l>=4.0*TK/CR»*2.0)»(TPC2>-TP(l)) 

DO 17 J=leICOL+l ^ 

FCICOB»$-»-6+J) = (F2(a)+RHO{J)»FlCJ)/iaO0.O) 

2 /CRX(J34RHO(J)»CP(a)) 

IFdCOND.EQ.O) F(ICQL46+6+J)=0.a 
CONTINOE 

RETURH 

':MO 
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55 

17 


16 


18 


22 

39 

23 

24 
26 


mm* 

OUTPU 

YC150) ,RH0(25) ,CP(25) , 


R,TK,U70,ICOWD, IYEJj,DR 


SUBROUTINE OUTPUT (XSOL , Y ) 

dimension - - 

I 

cSiSoN ICOL^c6Ll6f6^RX,UO,Tlif 
TYPE ♦.IT 

IFCIT.BQ.O) GO TO 36 
DO 17 Jaf,IC0l,+l 
sY(J) 

5¥C|CQLtl+J) 

= YaC0Lf2+2tJ} 

5YCIC0l,*3 + 3+j) 
=yClC0L*4+4+JJ 
aYCICOt^5+5+J) 
=YClCOL*6+6+j5-273. 
=Ry(4.J.IT)/RYC3,JiIT) 

RyCl.J.IT)/U10)»100.0 


n2C25J 


.(J)=UCJ)tRYC 6 :jilT> 
<Ji=utj)»CRYC7,J.IT) + 273.)*CPCJ) 


Ry(i,a,iT 
RY(2,j;iT 
RYCl.JilT 
Ry(4.j;iT 

RYce.aliT 
RY(7.J#IT 
RYC8,J.lTl 
RYC9.J4T) = £1, 

CONTINUE 
DO 16 Jsl.ICOL+l 

RHOC J)=1000, 0/(1 ,a065+0.0123*Y(J)+(RY( 7. J.TD +27 3.0-495. )♦ 

3 (0, 00035+0, 00007tY£j)3) 

CP(J)=0.6593jY(J)/Ui0 + Cl,-Y(vJ)/Ul0)l‘C486.1+0.337*£RYC?,J,IT) 

4 +273,0) )/S,00(i.£) 

COMTINUE 

DO 18 Jsl,IC0L+l 
U(J)=RHO(J)*RXCJ)>*'COL(J) 

UU(a)=U£J)»CP(J) 

UlCj)=U(J)fRY(l,J,IT) 

02CJ)=UCJ)^RYC2,J,IT) 

|l3(a)=U(0)»RYC3,J,lT) 

U4£j)=UCJ)l'RYC4lj,IT} 

U5(J)=UCJ)»i<Y(5;v7,IT) 

1161^*'- ■ 

U7i 

continue 
iPTsicor4+i 
IFAILi=0 

CALli DOIGAFCCQL.U.IPT.SUM.ER.IFAILI) 

DOlGAFCCOL.UU.IPT.SUMM.ERR.IFAlDl) 

OOlGAFkOLjOl , IPTIsUMI JerI i iFAIut ) 
DOlGAF(COL,U2.IPf,SUH2,ER2.lFAII.l) 
D01GAFCC0L,U3.IPT,SUM3.ER3,lFAll4l) 
D01GAF£C0L,U4,IPT.SUM4.ER4,IFAI01) 
D0iGAF£C0L;U5.1PT,SUM5.ER5.IFAILl) 
D0lGAF(C0L,06,IPT,SU«6,ER6fIFAIEl) 

^ ^ D01GAF(CQE.U7,IPT,SUM7,ER7,IFAIL1) 

IF(IFAIL.NE,03 GO TO 9§ 

UlAVsSUMl/SUM 

TYC10,IT)=CI.-U1AV/U10)»t00, 

TYC1UIT3=SUM4/SUM3 

TY<12.IT)=SUM6/SUM 

TYC13.IT)=SUM7/SUMM-273, 

TY(l4,If)-SUN5/SUM 
IFCXSOL.EQ.O,) GO TO 35 
WRITE€35|22) XSOU. 

FORMAT ( 1 1 X . * time ! ' .F 5 , 2 , 2X , ' ttr , % // ) 

NRITE(39,23) (RY(l,l,IT),I=l.lCnL+l) 

WRITE(39.39) CRY(2,I,IT),! = l,lCOli+l),CRYC3.I.IT).l = l,ICOI.+1 ) 
' ,£RYC4,I*IT)fI=lt 

i___ .t « £ ^ ... . _ . IG01.+ I) 

6 , / . 


CALL 

CALL 

CALL 

CALL 

CALL 

CALL 

CALL 

CALL 



A I i /AVriU,'?,/ 

2 9X. 'WATER :'.10F10,6 ) 

FORMAT (8X, ' MONOMER : ' , lOF 1 0 . 6) 


WRITE(39,a 

format C| IX 


-24) (RY(7.I.lT).Isl.lC0L+i) 

f »j5\;T4»i V * ix» Tempi .1 of 10 , 4 ) 

WRITEC39,26).CRY(8.I.lT).I = lj.IcgL+l> 
FORflAT£/,SX. 'CHAIN LTH.I *,iOFlO,3) 
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WRTTEC39,47) fY(10,IT),TyU3rlt),TYCtl,it)rTY(.14,lTj,TY(12,IT) 
F0RMAT€/,7X,'AV Coliv= Sf10.6,5X, 'AV TEMP= ',F30.4,/, 
t 2X,*AV CHAIR LTH= %F10.6, 4X , *AV DIMER= ',F10,6v/, 

2 6X,'AV RAfER= ',F10.6) 

IT=IT+i 

XSOL=X+Fl.OAT(IT)fH 
RETURN „ 

TYPE ♦^IFAIEi 
FND 

SUBROUTINE PEDERV(T,Y, PR) 

DIREMSION YC1503,PH(70,70) 

RETURN 

END 
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C 

c 

c 


17 

18 

19 

20 




11 


33 

C 


22 


23 


FINITE DIFFERENCE; METHOD H 


COMMON IC0L.c6D,fil6,RX,U0,TN,R,TK,U70,TC0Nr,TvFL,DR 
external FCft, OUTPUT Pe6erV > f f ' 

I ”°'lhltt 

3 Ul4fll20ju3O,U4O,U5O,U6O,O7O 

ICOL='NPTS-i 

IFCICOND.EQ.Ol NRlTEC39rl7) 

IFC|c0ND*EO.i5 WRlTE(39ll85 
IFC!cOND.EQ, 25 NRlTE(39ll9) 

FORHATC 15X, 50 C '*'),///, 30X, 'ISOTHERMAL REACTOR* ,///, 1 5X , 50C '♦' ) 

pRj|AT( 15X, 50 C**'),///,30X, 'ADIABATIC REACTUR ',///, 15 a, SOC *=t= * ) , 

pRMATClSX^S0C'#'||///^22X,'TUBUDAR REACTfjR WITH HFAT 

WRITE C 39, 20) 'RfNpTS^fjfojuloJCl^oJlJSO, 070,030,060 
F0RMATC5X. 'RA{)1US=' ,F6.2,5X, 'NO OF F.D. PUTNTSs ' , 13 , / , 

1 5X,'MOflOMER=*,F8.4,5X, '1st MOMENT^ ' ,F8 .4 , 5X, / , 

2 8X, 'PC1)=:',F8, 4, lOX, 'DIMERS', F8.4 si6x, 't£«P= * F10,4,/, 

3 2X,'0tft MOMENTS*, F8, 4, lOX, 'WATERS ,F8, 4,/, 7B('^'')) 


0.2) WRITE(39,21) 00, TM 

X,'iiEAT TRANSF COEFFs ' .FlO ,5 ,2X, 'KCAL/ (SO M 
',/,10X, 'COOLANT TEMps * ,F10.5 , // , 75 C ' ♦ ' ) ) 


TN®TW**273 

IFClCOH 0 ,feo. 2 ) WRITE(39 
FORMAT C/,5X ‘ ^ 

1 (DEG K) ■ 

TOLal.E-Oe 
IRELAB=2 
MPEDsO 
XsO.O 
0=0,5 
ITS0 

DR=1,0/IC0L „ 

DO 11 I=1,CIC0L+1 
CnLCI)=FL6AT(I-l)%DR 
DO 33 J=1,IC0L+1 

RXCJ)= 2 ,*Ci,-(COLC.jn** 2 .) 

IFCRXCjJ^EQ.O.O) Rxt0)=l.E-04 
IF(IVEL,eo.O) RXCJ)=1.0 
CONTINUE 

INITIAL CONDITION 
DO 22 I=1,(IC0L+1) 

YCI)=U10 
YCIC0L+1+I)=02« 
yaC0L*2+2 + I)=030 

¥aC0L*3 + 3 + I)=U40 
YaCOL»4+4 + I)sOSO 
Y(IC0L»5+5+I)=0C 
YCIC0L*6+6+I)=U1 
CONTINUE 
IFCIVEL.EQ.O) go to 23 
Y(ICaL+l)=0. 670647 
Y C ICOL ♦ 2 + 2 ) =0 . 00 0 2 3 4 
Y(ICOL*3+3)sO. 043585 
YClCOL*4+4)s8, 05944 
Y(IC0L»5+5)=0. 034956 
YaCOL*6+6)sO, 116415 
IFAILsCi 

NsCICOL+I)=»!7-2 

IWsiStN 

CALL D02EBF CX , XENO , N , Y , TOL , IRELAB ,FCN , MPF0 , 


) (HR) 



1 PEDERV, OUTPUT, WtlW.IPAIL) 

IFCIFAIL.NE.O 3 So TO 24 
STOP 

TTPE ♦,1FA114 


• TCOND, II/Eii,DR 


§ySHOyTlNE FCN 

|lSiNSl5NS(5|5|;|tl^O),AO(5),EO(S),FK(S,25i, 

COMMON ICoLCOI.,u!o,RX,u6,TM,R/fK,U70, 

DO 44 Jsl,ICOL«l 
TPCJ + U=YCjCOL*6+6+J} 

r ^ ^ 1 , > / n 

))»(TPCTCnL-13-4.0 

TYPI: ♦tTP(ICnL),TP(ieOL+J),TP(ICOL+2) 

DO 16 .J=l,ICOL+l 

RHf3CJ) = l0O0.O/(l jd065+0.0123*Y(J) + (TP(J)-495.)1‘ 

3 (0. 00035+0. 00007*YCJ))) 

CPCJ)=0.6593*YCa3)/U10 + C1,-YCJ)/U10)*(486.1+0.337»TP{J} 

4 )/10OO,O 
COfiTIMUE 

AOCl)=5.9874E+i3S 
A0(2)=1.8942S+10 
A0{3)=2.8558E+C)9 
AOC4>=8.5778E+li 
AOC5)=2.5701F+08 
E0C1)=1,988E+C4 
EOC2)=2.3271E+04 
E0(3>=2.284SS+04 
E0C4)=4.2E+04 
£0(5^2. |3E+04 
rlCl) = 1.9l8E+03 
S-5.9458E+03 
=-4.0438E+03 
S-9.6E+03 
="3.1691E+03 
“ 8846 


Oil Cl 
0HC2 
DHC3 
0HC4 
DHC5 
DS 


|) s >- 7 , Qg<iO 

2)=9.4374E-01 


03(3)3-6.9487 

DsJsJ^iiaSliE 


■01 


AC (1)34, 

AC(2)=1,- 

Ach)=1.6 




E+07 
,E+10 
7E+iO 


3307E+i2 

_.i|lE+09 

EC(i>3l,8§06E+04 
ECC2)s2.067E+04 
BC(33-2.dlO7E+04 
ECC4)=3.74E+04 
EC€5)32.04E+04 
UR31 987 
DO 11 lBl,5 

00 11 J31,IC0L+1 

FK(I,J)=AO(I)*EXPC-EOCl)/(UR»TPCJ)))+ACtI)*EXP 

1 (-EC€I>/CURnPCJ)3)*Y(ICOl.+2+2+J) 
EKClrJ^J'EXPCDSCD/UR-DHCD/dlR^TPCJ))) 
BKCl,J)sFKCl,J)/EKCl,a) 

CONTINUE 

00 12 J=1,ICQL+1 

P23Y(IC0l*+l+vJ) 

P33YCIC01.+1+J)^ „ ^ . 

APR=CY(ICOL4i2+2+a)-2,*Y€lCOL+l+a)) 
IPCAPR.OT.O.O) APR30.Q 
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C 

C 

c 

c 

c 

'Mi 

illl 

c 

12 


15 

C 

17 


IFCAPR.LT.O.O) P2=0.0 
IFCAPR.lt. 0.0) P3=0.0 


#li6o£J|+5||55|j?i*jf jf{35i$?ifoLl5+5+j)-8K(i;jHrcicoL+i+j) 

1 -FK€2,J)*YCICOL*2+2+J)»YflCOL*2+2+J)+ 

i „BK(2.a)J|aC0Lf5 + 5+a)*(y(TCDL*3 + 3+J)-Y(TC0L*2 + 2 + J)) 

3 +FKC4ja)*Y(IC0L*4 + 4+J)»Y(IC0L»5+5+J)- . 

4 iBKC4.J)tP2 

F(ICOL»2+2+J)=F(ICQL*2+2+U)/RXCa) 

fn6o£*K5;J5i«!i:j!H!JiH!icaE»5«;j!i6Kn;3!»i(lf0ttuj) 

I +FK(3,J)*y(J)»J(ICOI,»2t2+J)- 


f?ICOL*4+4+J)=-F'K(l,a5lY?|£o£^l+l+j)^Y?ICOL*§ + 5+vJ)+BK(i,J)*P2 

1 -FK(5,d)»YCICOLt|+4+vJ)iyCICOL*2 + 2+J) + 

2 8K(5,J)*APR 

r(XCOli*4+4+J)s:F(ICOL»4+4+J)/RX(a) 

?^^cgLj§+§+j5=-f'ftUppY(aHYct£6Ll5+§4j5^Kti;354ttic6L+i+a) 

1 +FKC2iJ)fyCICOL*2+2+J)»Y(ICOL»2+2+J)- 

2 BK(2tJ)=t'yClCOL»5 + 5+J)»(YaCOL»3+3+J)-y(TCOL*2+2+a)) 

3 -FK(4,5)*Y(IC0L*4 + 4+a)4YhC0L»5+5+0)+BK(4ra) 

2 ♦P2 

FCICOL*5t5+J)=FCICOL*5+5+J)/RX(J) 

cdSTiNSE’V’’’*'’'’'*’**’******’ 

DO 13 JJ=2,IC0L 
P2=YCIC0L+|+Jjr 
P3=Y C ICOL’*'! + JO 5 

^S?7O^lS°l*lt2t^^>r2.4Y(iC0L+l+Jj)) 

IFCAPR.LT.O.O) APR=0.0 
IFaPR.XT.O.O) P2*0.d 
IFCAPR,LT,0.0) P3=0.0 

FI C J J >=*gH ( 1 ) * { FK C 1 , J J ) ♦Y ( JJ 5 ^Y C ICOL^B+S-*- J J) «BK £ 1 , J J ) ^ 

1 -0H(25*(FK(2;jJ)»Y(ICOL*2+2+JJ)»Y(iCOL*2+2+JJ)-BK(2,JJ)* 


Y£ICOL»5+5+jJ)»(Y(ICOI.43 + 3+JJ)-YCTCOL»2+2+JJ))) 

-DH ( 3 ) » C FK ( 3 , J J ) ^^y C JJ ) *y ( IC0L*2+2+d J ) -aK( 3 , jj) ♦ 
YCICOL*2+2+JJ)+BK(3,JJ)=«'YCICOL+l+JJ}) 

-DHC4)»£FK(4, jJ)»y(IC0L*4+4+aj)*Y(lC0L*5t5+JJ)-BK(4.jvJ)*P?i 
-DHC5)’i<(FK(5,JJ)*YaCDL»4+4+aJ)*Y(lC0L*2t2+aj)-BKf5rJJ}*APR) 


CONTiritlE 

Do*i§*£=§nc6L 

F2CL)=CTK/R*»2,0)*C(TP(L+l)-TP(L-l))/2,a/DP 

1 /COI«CL) + (TP(L+l)-2,0»TPCL>+TP(L-l3)/£DR=^*2.O)) 
COHTINUE 

TYPE ♦»n(IC 0 L+t)#F 2 (IC 0 l-H) 

DO 17 Jsl-ICOL-l 

r C ICOL*6t® ♦ ^1= ^ P’2 i J-*- 1 f 1 f 1 ^ ^ 5 

2 /CRXCJtl3*RH0CJ+l)=*'CP(J+U) 

IFCICOHD.EQ.O) FCICOL*6+6+J)sy.O 
CONf IHUE . 

RETURN 

■END,.; 



75 


SySRQUIINE OUTPUT 

SUBnOUTIME 01JTPUTCXS0L,¥) 
dimension yC15C^),RH0C25>ICP(25), 

COgMON ICoLc6LM6Ux,UO,TW,H,TK,lJ70,TCnhP,TVEU,DR 

IFq|.|Q?0)^GO TO 36 

TP(J+l)sY(icOL*6+6+J) 

|pc||ot|!}=(uoi|w-(ipc2!6n)P*Rn»CTPCicoL-i)-4.o 
t ♦TPCIC0L3))/(U0+C3.0»TK)/(2.0»DR*P}) 

DO 17 Jal,lCOi4+l 
RYdril^ITjsTCJ) 

RY(2,J,IT)=Y(ICOL+UJ) 

RY(3,JrIT5«YqC0L^2+2tJJ 

Ry(4ra;iT)=Yhcoi,*3+3+j3 

RY(5W»IT)=YCIC0L*4+4+ji 

RY(6'jJlT| = YhC0L»5+5+j5 
RY(7,4,IT)=TP(J)-273. 

RY(8,J,IT)=Ry(4,J,IT5/RY(3,JtIT) 

RY(9|J^IT3=Cl,-Ry(l,O,IT)/U10)»lOO,O 

DO 16 J=l-ICOL+l 

RH0C0)=l0o0,0/(i,0065+0,0123*YCJ)+CRY(7,J,lT)+2?3,0-495.)» 

3 C0,00035+0,OOOD7*YCJ))) 

CP(a)s0.6543*Y(J)/U10t(i.-YtJ)/U10)»C486,l+0.337*CR¥C7,a/lT) 

4 +273.0))/1000,0 
COSTINUE 

DO |i asl,ICOli+l 

U(j5sRHQ(0)*RX(J)*C01,(J) 

UUCJ)SUCJ3»CP(J) 

Ul(a)=U(J3*RY(l,J,IT) 

02(J)=U(J3»RY(2,J,IT) 

U3(J)=U(J)»RY(3jjdT) 

U4Ca)=U(a34RY(4,0f,IT) 

U5{J)=U(J)»RYC5,J,IT) 

U6Cj5=U(j5*RYC6p,IT) 

U7C35sU(J)4(RK7,.5,IT) + 273.)»CP(J) 

CONTINUE 

iPTsICOL+l 

IFMH=0 

Ckhh D0lGAF(COL,U,IPT,SUIi|,ER,IFAXLl> 

CALL DOlGAFCCOL/lIU,IPT,SUMM,ERRrIFAILl) 

CALL DOlGAFCCOL,UdIPT,SUMl,ERl rlFAILl) 

CALL 001GAFCC0L,U24pT,SUM 2,ER2,IFAIL13 
CALL DOlGAFCCOL,U3,IPT,SUM3,ER3rIFAILn 


CALL D01GAF(C0L,U2,IPT,SUM2,ER2,IFA 
CALL DOlGAFCCOL,U3,IPT,SUM3,ER3rIFA_.... 
CALL D01GAFCC0L,U4,IPT,SUM4,ER4#IFAILI) 
CALL D01GAFCCOL,U5,IPT,5UM5,ER5rIFAILl) 
CALL DOIGAF CCOL , U6 , IPT , SUMS , ERS , IF AlLl ) 
CALL DOIGAF CCOL ,U7,IPTlsUM7,ER7rlFAIH) 
irCIFAIL.NE.O) GO TO 90 
UlAV=SnMi/SUM 

TY(10,IT) = C1.-U1AV/UIO)=I'100. 
TY<il,IT)=SUM4/SUM3 
TY(12#IT)=5UM6/SUM 
TY(13,IT3=SU«7/SUMM-273. 
fY(14,IT)=5UM5/SUM 
IFCXSOL.EO.O,) go to 35 
if RITE (39, 22) XSOL 

FORMATCnX, 'TIME: ',F5.2,2X, 'Hr, %// 
i#RITE(39,2i) (RYC 1 ,I,IT) , f =1 , iCOL+ 
ifRITEC39,39) (RYC24f IT) ,I = UICOL+l 

1 ,CRY€4,I,IT) .1=1/ 

2 ICOL+l5,tRYC§,I,I|),I=ltICOL+l)/ 
FORpTClOX/'PCl) J*/20FIO,6///5X/;0 
t 5X/*!st MO»ENT;%20F10.S,//9X, 'D 


/ (RYC 3 , 1 / IT) ,1=1 , ICOL+l ) 


ICOL+l) 

0,6/// 

/ f 
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23 

24 
26 
29 
47 


35 

99 

C 


2 9X,;WATER S%20F|0,6 
FOBMAI (SX.^MOMOMEH: %20Fl0.6) 

WRITE C3?f 24) CRY(7,I,IT),Iai,ICai, + l} 
FORMATCllx, 'TE mPj ', 20Flol4) 




iiPiPi;iSi;pSi: 




V= •,F10,6.5X,'»Y ie«p= 
= ‘,rtO.J,4x, AV DIMERS 


■fAU:'/, 


FORMAf(/#7X/'AV COftV* VFlO.e.SX, ^AV TEMP= 'rFlO. 

1 2X,'AV CHAIN I,|H5! %FIO,§,4x>*AV 0IMER= %F10, 

2 6X,*AV WATERS %F10r6) 

WRITE(39x27>,_ 

FORMAT(75C'.')) 

ITsiT+i 

XS0L=X+FI«0ATCIT)*H 

RETURN 

TYPE »,IFAIL1 
END 


SUBROUTINE PEDERVCT, Y ,PW) 
dimension Y(150) ,PWC70,70 
RETURN 
END 


) 



7 / 


C 

c 

c 

c 

c 

c 

c 

c 

c 


33 

59 

69 

37 


53 


20 

C 


G 

18 


f ♦ f f f f f ;^! 3(t * # 3(t * lj[ f ^ f jf, # ij! ^ Jj; Jj. ^ 

FIfllTE 0IFF; METHOD III 

TUBULAR REACTOR SIMULATION 
MAIN PROGRAM 

t^t*********’$;^*H(:^t$it:)lliilfi^^i^;^t***>¥t**i^*********$$**t 
MMENSION 017,50), U0UT(7,20),WRK (5000) ,X(5u},XOIIT(20} 

? ^|^i5)/J4C15) .¥5Ci5),Y6C15),Y7(15),RYC25r5?6,100) ,TYC25,100 V 

IxTiliimsIIJElf 

0pEN{UNl|s25»DEVlCE='0S{C',FILEs'flYL,0UT') 

OPEN (UNIT=f 5, DEVICES *DSK',FILE='’NYL, IN'} 

OPEN ( UNITS55, DEVICES 'DSK', FILES 'PLOT. OUT') 

NPDE57 

READ(15,*) ICONDrlVEL, 

1 UO,BB,TM, 

I B,NPTS,L,- 

LksSOOO 630,040, U50,U60,U70 

IFCICOND.EO.O) WRITEC25,33) 

IFfICOND.EO.l) >miTEC25,59> 

IFaCONp.EO-2) NRITEC25,69) 

FORMA|(15x, 50( ' =!='),///, Sox, 'ISOTHERMAL REACTOR ',///, 1 5X , 50 (' *') 
F0RMATCi5X,50( '* '),///, 30X, ' ADIABATIC REACTOR ',///, ISX , 50 C 1 , 

WRITE(25,37} B.NPTS, 010,040,020, 050,070, 030, 1160 
FORNATCSX. 'RADiys=',F6.2,5X, 'NO OF MESH POINTS* ', 13 , /, 
lX,:?jg60MER=' ,F8,4,§xl'l5t MOMENT*' eF8, 4, 5X,/ 


^K'.:fci)=:..Fa:4;!ox;'6iMER*'3 


moments *",F 8; 


M X -* ■ f r Q. ♦ *1 # j(. U.A. 

4,10X, 'WteR*%F8 


, T..:*’Pr * 

.4,/,7'^C' 


,F10. 

*')) 


4,/, 


3 2X,'0tL 

TN=TM-273, 

IFClCOND.EQ.p HRITEC25,53) 00, TN 

F0RHAT(/,5X, 'HEAT TRANSF COEFF* ' ,F10,5,2X, 'KCAL/f SQ M)CHR) 
1 (DEC K) ',/,10X, 'COOLANT TEMP= ' ,Fl0.5,//,75f ) 

M*! 

AsO. 

ACC=0.5E-06 

IMESH=4 

INITIAL CONDITION 
DO 20 I=1,NPTS 
U(1,I)*U1© 

U(2,I)=020 

U(3,I)*U30 

OC4,I)=U40 

0(5,I)=U50 

0(6,15=060 

UC7,I)=U70 

COHTINOE 

GO TO 18 

IF(fVEL.EQ-O) GO TO 18 
0(1, NPTS)=6. 670647 

S • 

0(4, NPTS)=8, 05944 
0(5, MPTS)=0, 034956 
0(6, »IPTS3*0, 116415 
U(7,NPTS3=TM 

TiPE *, ((0(I,J),a=l,NPTS},lsl,7,6) 

ts*o, . 

IND*0 

DO 60 IT=1,L 
T0Uf*FL0ATUT)/2. 

TX(IT)*TOUT 

IFAIL51 



78 


C 

65 

90 


13 

C 


r 

695 



C Ar.L DO 3PliF C gPDE , M , f DEF , BNDY , A ,B , TS , TOUT , U , fiPTS , IMESH , 
P0„13^J=1 .NPTS 

RHOslOOO^O^|1.0065|0,Ol23*UCl,J) + (UC?,J)-495.>>l=CO.f>OOJ5 + 

fa5/U10+(l.-UClrJ3/Ul0}*(486,l+0.3i7X‘nt7, j) ) 
TO 65 

WR|tEC24»*) RHO.CP 
IF(IVEI,,ftEa) g6 TO 90 

y3Cj3sYCvJ)*U(3 J) 

Y4Cj5»Y(a)*UC4,j) 

y5(jS=Y(J)*UC5'j> 

Y7(J}=Y(J)»CP*6(7,J) 

YY€J3-Y(J)#CP 

corivCJ)={i,-uci , j)/ui(j)=»'ioo, 

CMU(J)=:0(4,J)/0C3fU) 

TE.MPC7,J)=6(?,J)-i73. 

COWTI^UE* 

GO TO 695 

DO 63 J=l,(llPfS-l) 

HCJ)=XCJ+15-XCJJ 
OftV=UAV+0,S*HCj)*(YCJ)+YCJ+l) ) 

DUAV=OyAV+0.5!fHf J)*CYYCJ)+YYCJ + 1)) 


■'.r m w • * ^ w ^ ^ \ ’f-n* /f \ \ J 1 X \ hJ * X J 

U5=U5+e(J)*CYSCJ)+YS(J+l) j 
6=U6+C^5JH(an(Y6CJUY6(a + 13) 
!l7=U7+fj,5=<‘HCJ)*(Y7(J)+Y7(0 + l)) 

CONTIjIUE 


GO TO 64 
IFAliil^O 

CALL DOiGAFCX#Y,MPTSfUA¥,EP,IFAILl) 

CALL D01GAFCX,YYffiPTS,UUA¥,tRH,IFAILl ) 

CALL DOlGAF(X;yi;f^PTS,tJl,ERl.I#AlLl) 

CALL 001GAFCX,Y2,l«PTS,y2UP2aFAILn 

CALL p01GAFh,Y3 NPTS ul ErS IfAlLl) 

CALL DplGAFCX,Y4,f'IPTS,y4,EP4,IPAlLl) 

CALL p01GAF(X,Y5,NPT.5,n5,EP5,irAILl) 

CALL »01GAF(X,Y6,MPTS,U6,ER6,I^AIH ) 

CALL D0lGAFCX,Y7,MPTS,U7,ER7,I’^AILl J 

IFCIFIILI.NE.O) GO TO 99 

UlAVsUl/UAV 

U2AVSU2/UAV 

II3AVSU3/UAV 

U4A?=U4/UAV 

U5AVSU5/UAV 

U6AV=U6/UAV 

07AVSU7/0UAV 

U7AVC=U7A¥-273. 

COHAV=(l.-IJlA¥/UlO)»10O. 

X«U=U4AV/<J3A¥ 

|YC¥AaiABLE NC3,, IT) 

TYCI,,;) I I; 10=A¥ CONV ? ll=AV CHAIN 1 


10-A¥ CONV ? ll=AV CHAIN LENGTH 
12=A¥ WATER? 13=AV TEHPCC) 
i4=AV DIHER 


DO 28 Oai,NPTS . 
RY(7,J,IT)=TEMP(7fJ) 



28 


C 

C 

22 


39 

23 

24 
26 
29 
47 

27 


C 

691 

101 


103 

102 

104 

105 

107 

106 

108 
60 
19 

21 

693 

109 

C 

C 

C 

692 


Ty(l0,IT)=C0NAV 

TYCll,IT)sXMU 

TY{i2,n:)=u64y 

TY(!3,IT>=U7AyC 

TYa4.ITl=USAy 

WRIT|C25|22) toot 

FORMA|(llx, 'TIME: %F5.2r2X, 'Hr, %//) 
WRlTEi25,395 C0(2rI),Isl,RPTS) 


i__NPT5} ^**^fOF10^ * I £ Xrl 


F0RMAT(10x, *PC1 

1 5X,nst MOMEl 


ciK3,i),i=i»;'iPrs),(iic« 

rI)|I=l,WPTS) 

.,.6,/,5X, 'Oth MOMEMT; ',10F10.< 
2 :',10F10..,/, 

FORMAfjCSX, 'MONOMER: lOFlO, 6) 

WRI1’e125-24) CTEMP(7.I),TsLnptS) 
FORMAtCllx,'TEMP:M6Fl6.4r 
WRITE C25, 26 ),( CM0CI>,T=1 jNPTS) 

9^^ PLINTH,: %iOFio.3) 

WRITE(25,295 (COrjV(I) ,1 = 1 ,NPTS) 

F0RMAXC4I, 'COiiyERSIUR :',i0F10,3) 

«RlTE(2b/47) ^C0f||!kV,li7AV(:;XMU«U5A¥,U6AV 


FoWlAt C/;.7'P : Ay ’ c6nV='.' ' ^'f 1 b , § , §X t * A.V_ f EMP= : , Fig, | , / 


1 

2 


2X, 
6X* 


av Curjvr ' FlU,ft,5Xt'AV TEMP= 
'AV CHAIM L’fH= 'r#10.6,4X,*AV DIMER= 
'AV WATERS ',F10,6) 


rFli 


\u 


tYPE '♦jfOiJT^Cu"(J,Uf J=1»7,6) 

00 TO t>0 

# S(E if! f :(!!|! :|r *** 3^ f ](!# )(t f # 1^ f ^ ,tr % f ^ ^ JK j(t !(j jH sjj 

WRITE (25,101) TOUT 

FORMATCIIx/tIMf: %F5.2,?X,'Hr, ',//) 

WRITEC25,l62j (IJ(1 , 1) , Isl -NPTS) 

HRITEi25,lt)32 C U ( 2 , 1 ) , Isi ,fjpTS) , (U( 3 , 1) , I^l , WPTS) , ( U C ■ 
FORM AT C 1 6 X , f 1 j ' : ' i { 2^0 ! ^ / J 5 X ^ ' 6 tii ^ MOMENT : ' , 1 2F 1 0 . i 

FORMAfpX. 'MONOMER: ',i2F10. 6) 

WRITEC25j!o 4) CTEMP(^,I),Isi,NPTS3 
format (1 IX, 'TEMP: ',12F10, 4) 

WRITE(25,l653 ( Cm{i(I),I*1,NPTS) 

F0RMAT(/,5X,'CHAIN LTH.: ',12Fi0.3) 

WRITEC25,107) (CONyCD^IsllNPtS) 

F0RMATC4X. 'CONVERSION :’.!2F10,3 

/,U6AV 

. ^ . . ,-/AV TEMP= *,F!0.4//, 

„ LTH= ',f'lU,6,4X,%’.V RIMERS '^FlO.ii,/, 

2 6X,'AV WATERS %FtO,63 



WRITEI25^1082 


) 


FpBHA|(7U - ,, 

i^RlTE(25 19 j 

FORMATC/^, 'THE GRID POINTS APE(r IN P,);') 
lF(NPT5,Ep,12) GO TO 693 
WRITE C 25 , 2 1 ) C X C I ) , 1= 1 , NPTS > 
FORMATCiOX,6F8,33 
GO TO 692 

WRITEC25,109) CXCD , I=l ,NPTS) 
F0RMATC10X,12F8,3) 

* *6of PUT’Foi’PLoffiNS’ * •••••••••••*•••••••• 

WRIfE(5$,=l') CX(I),Isi,NPfl5,fTXtJ),Jsl,L) 


I) ,Isl , 


,I),I=1 

# / f 



BO 


t 


99 

100 


c 

c 

c 


34 

35 


WRITE(55,*) CCTYCt^J), 1 = 10,14), J=l,l4) 

TYPE ♦,IFAIL 
TYPE ♦,IFAIL,tOUT 
STQ-P ' 

END 

§yB8Qyj I E5^F 


lVEL,BB,TK,DD,B,IT,Ja,UlO 


495.)t(O,OO035+ 
6+Ci.-0X{l)/tI10)»{486,l+0.337f!IX(7)) 


OPEN C UN * A •>■ i j • uc. V 
READ (15,*) ICOMD 
ITER=lfER+l 
|YPE *,I|ER,T,X,»XC1),DUX(1) 

DO 35 I"l,6 

iFllVEL.EQ.C) GO TO 34 
C Cl 5=2.1 (1,-CX/B)*(X/B)) 

irCX.EO.g) CCl)=J.E-04 
O£I)=0,5E>2O 

CQNTINOE 

r"i;i3o27"icU?fffr-^ 

CPs0.^593*0X(i)/ai 
3 /iooo, 

? 1 1 C X/B) * CX/B 5 ) 
irCX.fcO.B) C£7)=l .E-04*RHO*CP 
IFOIVEL.EQ.D) C(7>=RH0*CP 
'G C 7 ) sTK 

A0fl)=5.9874£+05 
AOC2|sl,8942B+iO 
AD£3)=2.8558E+<)9 
A0(45=8,5778E+11 
AO(5)=2.S7eiE+fiH 
E0(J)=1.988E+04 
EC3(2)=2.3271E+04 
E0C3)=2.2845E+04 
E0C45=4.2E+D4 
E0(5)=2,i3E+04 
DHC1)=1.918E+03 
DH(2)=-5.9458E+03 
DH£33=-4.0438E+03 

DH(53=-3.1« 


i691E+03 
j 8 8 4 5 

OS£25=9,4p4EMn 
DS(3)=-S,9487 
DSC4)=-14,52 
DSC5)=S,8265E-01 

|^fi3=4:|0lg£;07‘ 

AC,.-. - 

ACC45=2 
AC (53=3 

EC£i)slI8805E+d4 
ECC23=3,007E+04 
EC(3)=2.0107E+04 
EcU) = 3.74E+04 
ECC5)=2,04E+04 
UR=1 .987 
DO 11 1=1,5 
FKCI)=A{K1)*EXP(' 
1 0X(7)))*UXC3) 


DSCl)=-7 
:2)=9. 


C23=1.2U4E+iO 
C35=1,6377E+|0 
.33e7E+!2 
#yiiEto9 


» ' S" • # • # «. .4K • » « # «> • • 


•E0(I)/(TJBtUX(?)3)+AC£I>*EXpC-ECCl)/(lJP* 


11 


EK{I>=EXFCDSCI3/UR-DH(I)/(IIR*0XC7))3 

BKa3=3FK£T)/EK(I) 

CONTIMOE 



81 


C 

C 

c 

c 

c 

c 


c 

c 

c 

C12 


55 

66 


C 

C 

C 


44 

C 

20 

55 


P2=iJXC2) 

Pi=UXl2)^ 

APR=(UXC33-2.»UXC2)) 

IFUPH.lt. O.S) APRsO.O 
p2=q.6 

IFCAPR.LT.O.O) P3=0.0 

r i * -.S j f X ? f .f f' A S ^ f( ### #,«,»,» a, 



3) 


E Je \ V j ^ -fv \. 4m .ft J — JT V, i ^ U 

iK ( 3 ) ♦P2 “FK C 5 ) tux ( 2 ) ♦UX ( 5 ) tBx C5 3 ♦P3 


f”5=i^^UK!!ll!yx555;Sif!I3*i;x!55iFK(25*uxnjJu;(5!; 

I BK(2)»UX(6S»(UX(4)-UX(3))+FK(4)»UXC5)»UX(6)- 

4 , ' oK.C 4 } ♦r'2 

ft455f^HH0ilH}JGrcl5:§Kti3^OH5Um3)Juxcniux(^3r 

1 BKC55*CpC3)-UXC2n+2,»FKtS2tUXC5)*IJX(3)-?,*ljf^C5 3« 

2 APR+2.»FK(4)*UX(65*UX?5)-2.»BK(4)*P2 

f?§Js“#KU)l»XC§3»IIX?b);BKU1l^2-#^!55*UXC5)*UXt3jr‘*’ 

i BKC5)*APH 

f cG}«^FKH5*5f|niQXfG3+BK(iUuXC2)+FK?25*G^n5*OXn3-* 

1 9K(2)*UX(6)»CUX(4)-UX(3) )-FK(43*UX(5)»UXC6)+BKC4) 

2 ♦P2 

h**:0Hn5^!^KCj 


1 -DHC 

2 -DHC 

3 -OHC 

4 «OHC 


iUSIili'jSIHj:! 


-fiKc!)iux(235 


SKC2)*UX(6)*(nxC4?-UX(3))) 
>+BKC3}*l 


♦ fFK 

!(FKC35*OXCi jJUXj|)-BKj|)*yX{3>+BKC3}*UXC2)) 


♦CFKU^ 

*CFKC5)*UX(5)#UXC3)”B«CC 


)*P23, 

)*APR) 


DO 12 1=1,7 

IF(OXCl).LT.O.O) UX(1)=8,8 
CONTXi'iUE 

IF(ICa?^D.EO,0) GO TO 55 
F(7)=F1*RHO/1000.0 
GO TO 66 
F(7)=0, 

RETURN 

END ' : 

♦ f f 4; 4: )>t ♦ f ♦ f f * f )(£ * f j(! :([ jj; * i(t f t ^ f ^ f. jjc f |s i; jj. 45 i| ^ 

SUBROUTINE BNDY 

subroutine BN0Y(NPD£rT,UX,I8ND,P,O,R) 

DIMENSION PCNPDE),QCNPDE),RCNPQEJ,UX(NPDE) 
0PENt0NIT=15,DEVTCf=^DSK*,FILE=*Mn4.Ifn 
READ(i5,») ICONl), rVEL,llO,TK,-3fi 
IFCICONO.EQ.U UO=O.Et0 


*,¥ \ ,3k V Mil'll ^ « J, / 

IFaSND.NE.O) GO TO 
LEFT HAND B,C, 

DO 44 I=1,NPDE 
P(I)=0. 

COMTINUE 

RETURN 

RIGHT HAND B.C,. 

CONTINUE 
DO 55 1=1,6 

QCI)=1. 

RCI)=0. 

CONTINUE 

0(75=-TK 


20 



!R(7)sO, 

RETURN 






e .3 




V, ■ i' *; 


OPTHOGONAL COfjLOCATIONs f«ETHOD I¥ 

f 1 ^ p » P3C , uo , TVi , R /f K , 07 0 , 1 cn. n . n^K 
EXTERMAG FCIsi, OUTPUT, PEOERV 
Ope - 1 ( ’JM I T= 1 5 , 0E \? TCE= ' 0SK % F I LEs ' ?JOC . I N M 
OPSM^MlfsjgJOEVIctsVolK'^ITlFs'MocIoUTM 
S88S58KII’^f#%X^CE=''DSK'.FIl[iEs'ERR.OOT') 
OPE?nUNIT=45jtOEVlCE=:'D,SK* ,FIL£=^COI,L.OO'/'' J 

READCis,*) i6oud,ivel,ico£./ 

1 UU,TK,Tl«, 

2 R,XEI«Or 

3 ■ Uia,U20,U3U,U4O,Opi^,UU0,U7U 

CAUL COLOOCCICOii,CQttA,B,s,T) 

TYPE, ■»# CCOt(I),i=l ..K^OL+n 
WRITECi 5 ,*HpUI,J 3 ,J=l,TCOL+tnT = l,TCOL 4 -n 
MRITEC45,»3 CuKI.a),vl=l .iCOB+l) ,T = l,TCnL + l ) 

TYPE 4.c«t(I) isLrcai.+!} 

IFCICOSD.EO.U) WRlfEC39,17) 

IFCICO’fD.EO.l) «RlTE(39ll8) 

IF(ICO110*EQ,2) WRfTEC39,!9) 

FOfWAT ( 15X, 50 €'’?'),///, 3dX, 'TSOTHERMAL REACTOR %///, 1 5X , 50 C "♦' ) - 
J //) 

FORMAT(15X,5OC/#»),///,30X, 'AUTASATIC REACTOR ',///, 15X , 50 C *«' V, 


i® v|'OR 4 Alci 5 X,|OC.^=t‘'),///,. 22 X,'TUOUL 
Ssill I mANSFER^,///.l5X,5uC'»'),//) 


AR REACTOR MITO OF AT 


4 , / , 


Tfjttf f|-»27'3' . ■ ' . . 

IF(ICOM0vIq.2) HRITEC3P,21) Uu,TO 

F0RMAT</#5X, 'meat TRANvSF COEFFs ' 

1 ( OEG S } ' , / , i OX , 'COOLANT TF. 1P= 

TOLai,E-06 

IRELAU=2 

MPED=U 

X=0,0 , 

Hat), 5 
IT=d 

00 33 J=leIcnL+l 
RXCJ)s2,*h.-CC0LtJ))»*2.) 

I F C RX C J 5 . EO . 0 . 0 ) RX C J 3 r 1 . E-04 
IFdVEL.EQ.O) RX(J)=1.C 
COWTIHUE 

INITIAL CQMOITIQII 
DO 22 I = l,(ICOL+n 
Y(I) = 0|(> 
y(ICOL+l+I)=U20 
YaCOL»2+2 + I)=U3C 

y(IC0L=*3 + 3 + I) = ll40 

YaCOL#4+4tI)=U50 

YCIC0L»5+S+I)=U60 

YClCOL*6+6+I)aU70/500,0 

COMTINUE 

IFClVEL,EQ,fj> GO TO 23 
YCICOL+1) =0.670647 
Y ( ICOL» 2 + 2 1=0 . C 0 0 234 
Y(IC0L*3 t 3)=0. 043585 
y(ICOL^4T4)=8.C5944 
yClCOL^5+5)=0. 034956 
YCIC0L»6+6)=0. 116415 
IFAILsO 

H3(IC0Lt1)*7-1 


j P 1 0 . 5 , ?X f ' KC AL / ( SQ M ) ( HR ) 
%F10.5,//,75f '*')) 


non 


26 

24 


10 

44 

m 


a 


XW=:18 + M 
Ti?J=TW/500,G 

call D02EBF(X,XEND,N,Y,T0L,IRELAB,FCN,MPF;D, 

1 , PEDERV, OUTPUT, W.TW.IFAIL) 

W«ITEC39,26) CC0LCi),I=l,k6L + l) 

FORMATC16x,10F8,3> 

IF(IFAIL,NE.O) So TO 24 
STOP 

type »,IFAIL 

“ND 

SUBROUTINE FCN 
SUBROUTINE FCNCT,Y,F) 

DIMENSION AC20.20),YC70),F(70),AnC5),EO(S), FK(5,1u), 

1 BKC5.l0),EKC5,10),DS(5),DHt5),ACC5),ECCS),CnL(lO).HHOf lOj 

2 ,CPCi0).BXClof,Bdo,20)'TP(l0j 
COMMON A,ICOL,COL#U10,B,RX,UO,T*J'!,R,TK,I 
SUMxO.0 

DO 10 tsl.ICOL 


I 


.U? 0 ,IGnND,TVFL 


SUMaSUM+AtiCOL+l .I3*Y{lCOL*6+6+l) 

tp(icol+i 5 =sC-tk/A* 5 um+ijo*tn) — 


IFCICOND.EQ.O) TP(lC0L+iysU70/500.0 
DO 44 Jsi.lCgL 
TP(J)aYCIC0L*6+6+a) 

DO 16 J=|.lCgL+l 

RHOCJ3 = l0OO,0/(i.0065 + 0.0i23=^YCJ) + (TP 

3 Jo.OOg 35 t 9 . 00 § 07 ?YC' 3 ? 3 ) 


/(UO+TK»A CTCOL+1 , ICOT.+l ) /K) 


CP£J)s6.6593*Y(J1/UJ0+(1, 
4 ♦500.0)/ 1000, 6 

CONTINUE^ 
AQa)=5.9874E+0S 
A0£2)5l.8942E+l‘} 
A0C3 )s 2.8558E+09 
AOC4)«8.5778E+li 
AO(5)=2.570IE+O8 
EO£l>«i.988E+04 
eO(2}=2.3271E+C4 
E0(3)=2,2845E+04 
EOC4>s4.2E+04 
E0(5)s2.l3E+04 
DHCl)=l,9i8E+03 
DHC2)=:-5.9458E+03 
DH(3)=-4.0438E+03 
DH(4)=-9.6E4-03 
OtiC5)=-3,1691E+03 
0S(l)=-7,8846 
DS(2)=9,4374E-ai 
0SC3)=-6,9487 
DSC4)=-14.52 
DSC5)=5.8255E-0l 


’(.T)*500,,.495,)* 

•Y ( J ) / U 1 0 ) ♦ C 48 6 , 1 +0 . 3 3 7 ♦TP 1 0 ) 


,«««.«# « m ♦' 9 ♦ f • 


, .• m m 9 m * 


5 JS s' f f ♦ f • t tl 3 • ? JS • • • • * 

AC(|)=4,|075E+07 
ACC2)-l,2U4E+10 
ACC3)=1.6377E+10 
ACC4)=2.3307E+12 
ACC5)=3.01lE+09 
ECC1)=1.8806E+04 
EC(2)=2.067E+04 
EC(33=2.0l07E+n4 
ECC4l=3C74E+04 
ECC5352.04E+04 
URal.98l. 

DO ll I=|,5_ 

Pk ci 1 JO^wl I J ^lipC-EOC 1 )/<UR?TP ( J) *500 . ) ) + AC ( I ) *EXP 

BK£lfJ 5 *FK£I,a}/EK(l,a) 

COMfIMUE 

DO 12 J=1 ,IC0L+1 



85 


G 

12 


C 

14 

13 


P2=YCIC0f.+l+J) 

P3=:YCIC0L+1+J) 

APR=a(ICOL=^2 + 2+J)-2.*Y(ICOL+l+J)) 

IFCAPR.LT.O.oi ApR=0,0 
IFCAPR.LT.O.O) P2=(>.0 
IFCAPR.LT.O.O) P3=0.0 

FCj5=-FKa, J)fYCJ)*y acOL^5t5+j)+BKCl f C ICUL+l+JJ 
1 -FKf3,J)*Y(J)*^CICOLf2+2+J) 

FCa)=FCJ)/RX<J3 

FtiaaaiiaS=F^u;j5*!?jni?ia6£i5;5;j)-BKh%\b**YriaoL;^ 

i -2,*FK(2,a3»Y(ICOLtl+J3*Y(ICOLt2+2+J3 
I +2^||K|2|J|||lICQlj5|5+a)nY(ICOLt2 + 2+a3-Y(IcnL+l+J)3 

4 "*^BKcSta3Jp2-FKC5tjJJYilCOL + l+J)»Y(TCOL*4 + 4+J)tBK(5,J)*Pi 

FClCOl.+ l+a3sFCICOL+l+J5/PX(J3 

t§+5 I -B^ Cl, J ) *Y ( iCOlt 1 + J J 

2'i'2+J ) f , 

W.f\ V' A ^ Jt L 4 w T ;2 TU / ^ I, V lCOL*3+i4-J3-Y(ICnL4=2f24-J3) 

3 +FKC4.a3*YClCOL»4+43-J)»yCICnL*5 + 5+J3" 

4 BKC4.J3*P2 

FCICOL»2+2+j3=F(lCOL*2+2+J5/PXCJ3 

^tt354y?icoCi5+5+j)-bKH 




3 


U«\ 1. • w ’ V J. V A VWJU • fcT*XW / * C ICOL 4 1 + vJ 3 3 + 2 « "^^FK ( 5 » J 3 ^ 

YClCOL*4 + 4+33»YCICOL*2+2+il3-2.KBK(5£J3* 

.. APR+2.*FKC4*a3*YClC0L*5+5+J)»Y(IG6L*4 + 4+J)-2.*BKC4,J3*P2 
(IC0L»3+3+J3=F(IC0L*3+3+J3/RXCJ3 

#?i56£*IUW5=i?-KU:33*Yaa6£iU4+j3HUeQf.l§l§;jHMi3,jj*^^^^^ 

1 -FK(5,J3*YCICOl.*4+4+J)4Y(ICOIi*2+2+J3 + 

2 BKC5-J)4APR ■ 

F(ICOl.*4+4+J3=F(lCOL*4+4+J)/RX(J3 

Fcic6£*5+5+j)=-FKH #55 *Y!j344!i2of, 45+5+33 +PK( 1 135 Jycicol+i+u) 

1 +FK(2jra3*Y(ICOL42 + 2+J34Y(lCOL*2 + 2+33- 

2 BK(2,J)4YCICOL45+5+J34(YClCOL43+3+J3»yCTCOL42+2+J3) 

3 -FK(4#J3*Y(ICOL44+4+J3+Y(ICOL45+5+J)+BKf4#J) 

2 4P2 

F(IC0L45+5+J3«F(IC0L45+5+J3/RXC3) 

continue 

DO 13 Jj=l,ICOL 

P2sY(ICOLtl+JJ) 

P3=Y(ICOL+ltJJ3 

APR=(YClCOL»2+2+JJ3-2,4Y(ICOL+l+Jj33 
IFCAPH.LT.O.Oj APR=0,0 
IF(APR.Lt.0.03 P2=s0,0 
IF(APR,I#T,0.03 P3=:0,0 

Fl=-DHh3*CFKCl,JJ3 4Y{JJ)4YCICOL*5f5+JJ)-tiKCl#.Iu3 4YCICOL+l+JJU 
1 -DHC234CFK(2#a33*Y(ICOL42+2+333*lf(ICOI.42+2+JJ3-BKC2,JJl4 

‘ yCIC01.4S+5+j33*CYCIC01*3+3+JJ3"Y{IC0L*2+2+J33 3 3 

-DH(3>*(FKC3,JJ3*Yt333»Y(ICOlj»2+2+3J3-HKC3,,3J3» 





9 • • il * f) « « « # f H 9 , 9 . 9 ' ^ 


2 

3 

4 

5 

6 

SUM2=3l0**‘’'’ 

*2 ! f ill Jm 2?I (]a , ICOL+13 *TP 

FClCOL*6+6+3J) = (F2+RHO(a334Fi/lOOO.0/500,03/CRXCiJd> 

2 4RHOCJJ34CP(aj33 

IFCICOND.EQ.O) F€ICOL46+6+Ja)sO,0 

CONTINUE 

RETURN 


4#4J5*P2' 

5#J3)4iiPWi 
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C 

c 

c 


55 

17 


25 


16 


18 


19 


22 


|< f ♦ |! f f f ♦ ^ ^ 4i * 4; f f :([♦ j|! t f ♦# t f ♦ 

SUBROUTINE OUTPUT _ 

t*t*^*****il^***tiifi$^*ittt^i>t^$}^***^'^*tif^**^^**t^* *$t*4*1t'*** 

SUBROUTINE OUTPUT (XS01,Y) , , ' , , 

DI'4ENSION flC20,20),B(20,k),Y(70),WTa03,RHyClo;,eP(ia)/ 

COmSioN A?I?OL;cb£^6l6^i,PX,UO.TIl,R,TK,U70,ICnND,TVEL 

TTPF ♦ - IT 

IFCiT.iQ-O) GO TO 36 
DO IT J=1,ICQL+1 

:YCJ) 

:YCIG0L+1+J) 
lsy(ICOL*2+2+J) 

5YCICOL*3+3+j 5 
|i!y(ICOL=»'4+4+J) 
isY(ICOL* 5+5+J3 
l=YCIC0L»6+6+J)*500. 0-273. 
(=RyC4,J,IT)/RY(3,yj-IT3 
iaCl,-^YU,0/TT)/U10)=('100.0 


RYC1,J,IT 

RYCa.jIlT 

RYprAJrlT 

RYUrJfl'J 

RYCs'vJilT 

RYC6,J>1T 

RyC7,Ji|T 

RY(8,a,|| 

RYC9,J,If 
CONTINUE 
SUMlaO.O 
DO 25 K= 
SUMl=SUM 


. I 


OL 

ICOL+1 


o u iM 1 — Q u ra i T ft v ji u jj T * » K 3 ^ Y C J C 0 L 6 + 6 K ) 
fllAlT=C-TK/R*S0Ml+U0*TW3/(U0’i-TK*A(IC0L+l , ITGL+l )/R) 
IF(|COND.E0.03 TLASTsUTO/SOO.O 
RY<7 , ICOL+l ,IT3-T1AST»500, 0-273 .0 

RHO( J3^lisi?i/tl. 0065 + 0. 0l23*Y(J) + CHy (7/0^13+273.0-4 

aPCJ3so!65l3»y(33?U10+ci!-YCJ3/U103=«‘C486.1 + U.3B7 + CRYt 

4+273.033/1000.0 

SUm|=o!!I;SUM2=:G,?SUH3sO,jSUM4=0. ?sUm 5=0. ?SUM6s:0.J SUM7 = 0, 
SUM=0,;S»MW=0, 

DO 18 J=l,lC06+l 

U(J3=RHOCJ)*RX(J3 

UU(J3=UCJ3*CP(J) 

Ul{J3=U(J3+RYCl,vJrIT3 

U2(J3=sU(J3*Ry(2,JrIT) 

U3CJ3=UCJ)+RYC3,J,IT) 

U4(J3=UCaJ)*RYC 4,J,IT3 
U5C03=U(Jj+RyC5, JrIT3 

e!J:J?=Sb"j:?K*I!l3+273.3+CPCU7 

CONTINUE 
DO 19 J=1 ,ICOL+l 
SUM=:SUM+HT(J3’''U(J) 

SUMM=5UMM + WTCvjj+UUCJ 
SUMl = SUMl+WT(J3i!JlCJ, 

SUM2=SUH2+WT(U3*U2CJ, 

SUM3=SUM3+HTQ)*U3p 
SUM4=SUH4+NT(J3*U4p, 

SUM5=SUM5+WT(J3*U5CJ, 

SUM6»SUM6+WTCJ3*U6Cq. 

SUM7-SUM7+NT(J)*U7(a) 

CONTINUE, , ; 

UUVsSUMl/SUM 

TY < 10 , 1 T 3 5 ( I , -U lA y / U 1 0 3 ♦ 1 0 0 , 
fYUl»TT3=SUM4/SUM3 

tYU3llT3=SUM7/SUMN-273. 

TYC14,IT3=SUY5/SUH 
IFCXSOL.EQ.G.) GO TO 35 
WRITE(3ff2§3 XSOL. „ v 

FORMATCiix, 'TIME: '/FS.a^aXr 'Hr, ",//). 


K « A 1 ’ i i I A I i JL r- ; ^ # r j * ,4 r 4 4 f n r . . 

WRITE(39,233 (RYP4^ilt34l7Ul£P^tP rovf ^ t r+l i-l fCnU+iO 
WRITE (39 1 39) CRYC2rI,lT), ial»ICOt+l 3 # CBYC3f I # IT) * I- 1 # iruii+i ) 



8 7 


39 

23 
2 4 
26 
29 
47 


27 

36 

35 


C 

33 

C 

C 

C 

c 


c 

c 


I fcOL+ljJ jRYl^!!lifT),I=l,ICOL + l),(Ry( 

FflRMlT f I 0 f f 'V' . 1 Ac* 4 A & / M 


F 0 RMp 510 Xv*P(l) : * /10F10.5,/,5X, '0th 
J 5X,*lst MOMKNT: lOFiO,6!/.9x # 'dimer 
2 9X^n«T|R^j',l0F{0,§_)* ' ^ 


6,X,iT) ,1 = 1 ,lCuT.+ t j 

^ * * ... 1 A L'* i f'j y 


moment: ', iOFlo 


iOFlO.C),/, 


. tJ 


■ 10 . 


6) 

l,ICOL+l) 


FOHMATC8X, 'MONOMER: ',10! 

WRITEC39,24) CR¥t7,lI!T 
FORMATClIXj 'TEMP: ',10F10,4) 

IPIf /:|x^^Sj(P'JT;S!l-!:} 6 SS 5 i’ 

F0RflAl(4ii, 

WRITEC39, 47) |yC10,IT),TY(13,lT),TY(ll ,IT),T7(14,IT} ,TyC12,IT» 
FORMATC/r7X,'AV CONVs ' ,F10.6, 5X , ' A V TEMPs ' , FIO 1 4 , / , 

1 2X,'AV CHAIN L|Hs',E| 0.6, lx, T aV DIMERS ',F10. 6,/, 

2 6X,'A¥ WATii:R= ',F10,6) 

NRITEC39,27) 

F0RMATC75('.')) 

IT=IT+1 

XSOLsX+FLOAT(IT3fH 

RETURN 


RE TURN 
END 

SUSROpiNE COLLQC€N,X,A,B,W) 

IMPLICIT REAL >l‘8 CA-H,0-z5 

dimension «3inv(|o#10),zC10),c(lO,iO),dC10,10) 
dimension a(20,20) ,6(20,20) ,q(20,20) ,x (20) ,w{20j 
DIMENSION DlFl(lO),0IF2(lO),DIF3ClO),RnnT({o) 
COMMON/ONE/ X1C2<>) 

N=2 ■ 

NTsii+l 

BETAsO.O 

aa=2, 

CALL JCOBI C 10 , H , 0 , 1 , 1 . ,BETA , D IFl ,DIF2 ,DIF3 , ROOT) 

ROOTCWT)=1.0 

DO 33 Kssl ,N 

Xl(K)=(SQRTCROOTCK))) 

TYPE ♦,Xi(K3 
CONTINUE 

call collCa,b,q,x,w,20,N,aa) 

TYPE ♦,(X(i5,I=1,NT) 

TYPE *,C(A(r,J),J=l,NT),ls!l,NT) 

TYPE ♦,C(B(I,J).Jsl-NT),l5l,NT) 

TYPE ♦,(W(I),Is{,MT) 

RETURN 
STOP 
END 

subroutine JCOBI (ND,N,NO.Ni,AL,BE, 

IMPLICIT REAL »8 CA-H,6-z5 

dimension DIF1(1O),0IF2(1O) ,DlF3clO),ROOT{iO) 


DIFl ,D1F2,DIF3,R00T) 




ROOTS AND DERIVATIVES OF JRACOBX POLYNOMIALS 
; MACHINE ACCURACY 16 D; 


FIRST EVALUATION OF COEFFICIENT IN RECURSION FURn:j:AS 
RECURSION COEFFICIENTS ARE STORED IN DIFi AND JTF? 


AB=AL+BE 

ADsBE-AL 

dIf I f tfs € ( AD/ C AB+2 . ) ) + 1 . ) /2, 

lF(l!ii?*2) GO TO 15 
DO 10 Is2,N 

Z=ABt2*Zl ■ 



BH 


it 


10 

c 

c 

c 

15 

25 


DIFtCIJs(((AB*AD)/(Z*(Z+ 2 *)J)+l.)/ 2 . 

IF Cl .WE. 2) <3Q TO U 
^DlFaCIJsCCAB+AP+ZD/CCZ+D^Z+Z)) 

GO TO 10 

z=z*z 

lf=Zl*{ABfZl) 

YsY^CAP+Y) 

0IF2CI3=Y/C(Z«13KZ} 

CONTINUE 

ROOT DETERMINAATION BY NEWTON METHOD WITH SUPPRESvSiO.M OF PRFViMrjs 


X=!0. 
DO 20 
XDSO, 
XN«1. 


I«1,N 


30 


22 

21 


20 


C 

C 

C 


57 


30 




)+X)*XN-DIF2(a3»p 

■ IF2CJ)*XD1+XN 


XD1=0, 

XfllaO. 

XP1b(*JiFK j)+X)#XNl-Dl 
XD-XN 
XDlaXNl 
XN=XP 
XNlsXPl 
CONTINUE 
ZC=*1. 

Z=XN/XNl 
IF (I ,EQ. 

DO 22 Js2,I 


1) GO TO 21 


ZCs2C-Z/CX-R00T(J-l)) 

ZaZ/ZC 

XaX-2 , ^ 

IF (DABSCZ>.GT,l.D-0g3. GO TO 25 

!rCA8SCZ).GT,l.E-06) GO TO 25 

ROOTCpsX 

XaX+,0001 

CONTINUE 

DIF3CI3=0. 

RETURN 

END 

SUBROUTINE COLLCA,fl,0,X,W,ND,N,AA) 

IMPLICIT REAL »8 (A-H,0-z5 ' , 

dimension ACND,ND),BCNDeND) VQCND.ND) ,XfNn)fVi(ND) 
dimension OINVC10,103,7C10),CaO,10)»D(10,10) 
COMMON/ONE/X1C20) 

NlaN+1 


THIS subroutine COMPUTES THE 
USING SYMMETRIC POLYNOMIALS 
DO 57 1=1, N 
XCI)=X1CI) 

CONTINUE . 

XCNDsi.O 
D0_30 1=1, N1 

zfifal./(2.*AI+AA-2.) 

continue 

DO 35 1=1, Nl 


MATRICES FOP orthogonal COLLUCATIOM 


ii*!!*;! 


;k 


0(I^J|=XCl5*»C2»J-2) 


J) 

0 40 

CA=2 . » J- 


KI,J3 
il : 


35 



89 


DA=C2.»J-2,)*C2.*JtAA-4.) 
DO 40 l=t,Nl 


40 


45 


50 

C 

12 

13 

C 

C 

c 

c 

55 


5 

10 


15 

20 


5 

10 

15 

20 


^ SEV X',— i. f i 

CCI,J)=CA*X{I)**(2*J-3) 
DCl,J) 3 DA*X(I 3 *=«'b»a- 4 ) 
COffTIMUE 


nn^kn^Tl? CpNV,Nl,i 03 
DO 50 1st, Ml 
WCI)s0,0 

no Ni 

||| 5 fg 5 ?IE 5 i^JxNVC.,X, 

contimue 

FORMATC'Bia VALUES ARE ') 


WRITE (60, 123 

FORMAT C9Et6,5/3 

RETURN 

END 

SUBROUTINE IN VRf A, N, MI) 

IMPLICIT REAL *8 lA-H.O-Z) 

D I ME N SIGNAC f i I , U 1 3 ,8 (203 ,0(4003 

DO 5 «3s I , N 

INDsN»Ca-13 

DO 5 Isl,N 

CCIN0+I3=A(1,U3 

CONTTNUR 

CALL irflERTCM,l,C,B ,13 
DO 20 Jsl,N 

00 10 Isl,N 
B(I3=0. 

BCJ3=1, 

CALL IMVSW(N,1 ,C,B,13 
DO 15 Isl,N 
A(I,J3 sBCI) 

COIJTIMUE 

CONTIMDE 

RETURN 

END 

SUBROUTINE INyERTCN,NE,A,B,ITYPE3 
IMPLICIT REAL *8 (A-R,0-b 

D IMEN SION A C 1 000 , 1 3 ,a C N 3 , A I ( 1 000 3 , B 1 C 1 000 3 ,G I C 1 000 3 
GO TO (5 , 10 , 15 3 , iTYPE 
CALL DECOMP (N, A) 

RETURN 

Mps(M-13/NE+l 
PAUSE 14 

CALL LUDEC0CNF,NE,A3 

RETURN 

DO 20 KsJ,N 

Al(K)sACK,l3 

B|{K|sAtK,23 

ClCk3sAiK,33 

CONTINUE 

CALL INVTRlCN,Al,8i,C13 
.DO 25 Kst-M' 

ACK,13=a|CK3 

ACK,235B1CK5 

ACK,33=C1(K3 
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CONTIf^UE 
RETURN 
ENTRY INVSW 

GO TO ( 30 , 35 , 40 ) fITYPE 
CALL S 0 LV£(N;A,B) 

RETURN 

CALL FAS(NP,NE,N,A,a) 
RETURN 

CALL 5 WEEPCN,A 1 ,B 1 ,C 1 ,B) 

RETURN 

END 


SUBROUTINE SING(I) 

IMPLICIT REAL fe CA-H, 0 -Z) 

GO TO ( 5 , 15 ), I 
piTE( 60 ,l 0 ) 

FORMATt 'MATRIX WITH ZERO ROW DECOMPOSE') 

RETURN 

formIt t^llS^ULAR MATRIX IN DECOMPOSE. ZERO IN SOLVE 
RETURN 
END 


subroutine DEC 0 MP(N,A), 
implicit real *8 CA-H, 0 -Z) 
dimension A(N,N) 
common /DENSE/ IPS( 201 ) ,SC( 201 ) 


IF (N .EQ. 1 ) RETURN 
DO 25 1 = 1 , N 
IPS(I)=I 
ROWfJRM= 0.0 

DO 10 0 = 1 , fi , _ 

IF (ROWWp-ABS(A(I,J))) 5 , 10,10 
ROWNRM sABS(A(I,J)) 

continue 

IF (ROWNRM) 15 , 20,15 
SC(I)=i,/ROWNR« 

GO TO 25 

CALL SING(l) 

SC(I)= 0 . 

CONTINUE 

MH 1 =N -1 

DO 65 K= 1 ,NM 1 

BIG= 0 . 

DO 35 I=K,N 
IP=IPS(I) 

S 1 ZE=ABS(A{IP,K)»SC(IP)) 

IF (SIZE-BIG) 35 , 35,30 

BIG=SIZE 

I 0 XPIV=I 

continue 

IF (BIG)^ 45 , 40,45 
CALL SING( 2 ) 

GO TO 65 

IF (IDXPIV-K) 50 , 55,50 
J=IPS(K) 

IPS(K)=IPS(IDXPIV) 

IPS(T 0 XP 1 ¥)=J 

KP=IPS(K) 

PIV 0 T=A(KP,K) 

KPlsK+t 

DO 60 I=KP 1 ,N 

IP=IPS(I) 

EM=-A(IP,K)/P 1 V 0 T 
A(IP.K)=-EM 
DO 60 J=KP 1 ,N 



A(IP,J)=aCIP,J3 +EM*ft(KP,J) 

GOWTIMUE 

continue 

KPsTPSCN) 

CONTINUE 

RETURN 

ENO 


SUBROUTINE SOUVECN,A,Bl 
IMPLICIT REAL *8 (A-H,0-Z3 
DIMENSION A(N,N),BCN) 

COMMON /DENSE/ IPS (20 13 ,SC(201 ) 


GO 

»13 


TO 5 


IF CM ,GT. 1) 
BC13=8C13/AC1 
RETURN 
CnUTINUE 
NPi=N+l 
IP=1PSC1) 

SCC1)=:BCIP) 
no 15 I=2^N 
IP=IPSCI) 

IM1=T-1 

SUMsO, 

go 10 j=i,iMi 

SUM=SUM+A(IP,J)»SC(J) 

CONTINUE 

SG(I)sBClP3-SUM 

CONTINUE 

IPsIPSCN) 

SC(fO=SCCM3/AClP,N3 

DO 25 IBACK=2,ri 

I=NP1"IBACK 

IP=IPSCI3 

IPlsI+1 

SUHsO. 

DO 20 J=IP1,N 

SUM=SUM+ACIP,J)*SC(J) 

CONTINUE 

SC(I) = (5CCI)-SI!M)/A(IP,I) 
CONTI HUE 
00 30 1 = 1 .N 
B(I)=SCCI) 

RETURN 

END 


SUBROUTINE FASCNP,NE.NT,A,B3 
IMPLICIT REAL *8 tA-B,0-z5 
DIMENSION A (NP , NP , NE> ,B (NT) 

NPlsflP ^ 

DO 10 L=1,NE 
DO iO I=2;NPi 
I2=(L-13»(NP-1)+I 
S=0 , 

11=1-1 
DO 5 J=t#Ii 
J2=I2-ItU 

S=S+A€liJ,L3*BC023 

CONTINUE 

B(I2)=eCI2)-S 

CONTINUE 
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B(NT)=BCNT)/ACNP,NP,NE3 

f<l=HP-l'' 

DO 25 K=1,N1 
I=N1+1-K 

I2 = (L**l>»CNP-l)tI 

M=I + 1 

M 2 sNi+l 

SsO, 

DO 20 J=M,N2 
J 2 e(l 4 -i)»CNP- 13 +J 
S=S+A(IjrJ,L)»B(J23 
B{I2)=CBCf2)-S)/ACl,I,L) 
GONTINUE 
RETURN 
END • 


SUBROUTINE LUDEeO(NP>NE,ft3 
MPLICIT real *8 (A-H,0-Z) 
IMEMSION ACNP,Np,NE3 


M 1=1! P-1 
no 10 L=1,NE 


no 5 K- 


: » n 

= 1 ,!il 


K 1 =K +1 
DO 5 I=K1,NP 
S=ACI,K,L)/A(K,K,L) 

A(I,K,L)=S 

DO 5,a=Kl ,NP ^ , 

ACI,J,L)=A(1,J,L)-S*A(K,J,I,) 

CONTINUE^ , 

IF CL,E0,NE> RETURN 
AC1,1,L+1J=A(NP,NP,L) 

CONTINUE 

END 


SOBHOUTINE IN VTR I CN, A , B ,C 
IMPLICIT REAL *8 CA-H,0-2 
DIMENSION ACN) ,B(N) ,C(N3 


DO 5 L=2rM 

S=ACL)/B(L-1) 

BCL)=BCL)-.S*C(L-1) 

ACL)=S 

RETURN 

END 


SUBROUTINE SWEEP C N , A^B.C , D ) 
IMPLICIT REAL *8 CA-HjO-Z3 ^ 
DIMENSION A (N3,B(N3»C(N),DCN) 

DO 5 L=2|N 

D(L)=DCl5-ACL3*D(L-1) 

DCN)sDCN)/8(N) 

00 10 L=2,N 


K=N-L +1 

D(K)sCD- 

RETURN 

END 


CK)*C€K>*D(K+1))/B(K3 



